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ABSTRACT 


Relativistic  neutral  particle  beams  (  Nl’Bs  )  off'c  r  the 
potential  for  highly  lethal  space-based  weapon  systems  and 
decoy  discriminators.  These  uses  of  NPBs  have  motivated 
(■nt.cri i  vi:  theoretical  and  experimental  studies  of  nega  I  i  ve 
ion  sources.  Negative  ion  sources  are  divided  into  two 
categories:  surface  sources  -  such  as  the  Penning  discharge, 
and  volume  sources  -  such  as  the  magnetic  multicusp 
discharge.  The  negative  ion  beam  required  as  input  to  a  NPB 
system  should  be  of  high  current  and  low  emittance.  Surface 
sources  produce  a  high  current,  high  emittance  negative  ion 
beam;  whereas,  volume  sources  produce  a  low  current,  low 
ein  it  l.nnce  beam.  Due  to  the  use  of  reactive  elements,  such  as 
cesium,  in  surface  sources,  they  have  a  much  shorter  useful 
lifetime  than  volume  sources.  While  current,  limitations  may 
preclude  the  use  of  volume  sources  in  NPB  weapon  systems, 
they  are  excellent,  candidates  for  use  in  NPB  discrimination 
systems.  Thus,  a  typical  volume  source  was  chosen  as  the 
topic  e  f  l,h  is  thesis. 

This  study  considers  negative  ion  production  in  a 
magnet  ie  multicusp  ion  source.  A  self-consistent  model  based 
on  Ma-.welJian  moments  of  the  Boltzmann  equation  is  presented. 
The  departure  of  the  electron  energy  distribution  from  a 
Maxwe| I  j an  is  also  discussed  and  incorporated.  The  results 
of  this  mode]  for  discharge  currents  ranging  from  1-100  A, 
for  pressures  ranging  from  1-100  inTorr,  and  for  voltages 


tang  ini'  from  20-100  V  are  compared  to  previous  studies  and 


v 


used  to  formulate  scaling  laws.  The  scaling  laws  predicted 
for  the  negative  ion  density  are  seen  to  be  supported  by 
previous  experimental  observations.  Finally,  the  scaling 
laws  are  used  to  predict  the  optimum  operating  regime  for  a 
magnetic  multicusp  ion  Bource. 
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CHAPTER  I  -  INTRODUCTION 


The  use  of  neutral  particle  beams  in  heating  of  fusion 


plasmas  and  in  weapons  research  has  motivated  considerable 


interest  in  the  study  of  sources  that  produce  beams  of 


negative  ions  with  high  brightness  and  low  emmitance.  One 


such  source  of  interest  is  the  magnetic  multicusp  ion  source. 


Many  experimental  and  theoretical  studies  have  been  conducted 


on  this  type  of  source  in  order  to  determine  the  processes 


involved  in  the  volume  production  of  negative  hydrogen  ions 


and  in  order  to  optimize  the  ion  production.  Before  these 


studies  are  discussed,  a  typical  magnetic  multicusp  ion 


source  will  be  described. 


The  Magnetic  Multicusp  Ion  Source 


The  magnetic  multicusp  ion  source  (MMIS)  was  originally 


developed  to  improve  the  containment  of  a  collisionless 


quiescent  plasma  [1].  A  source  of  this  type  consists  of  a 


;tal  container  with  some  wire  filaments  on  one  wall  (see 


Fig.  1).  The  container  is  filled  with  low  pressure  H2  gas 


(1-40  mTorr,  typically)  and  a  current  is  passed  through  the 


filaments  causing  thermionic  emission  of  electrons.  The 


filaments  (cathode)  are  typically  biased  negative  with 


respect  to  the  chamber  walls  (anode)  in  order  to  accelerate 


the  emitted  electrons  through  the  hydrogen  gas,  creating  a 


p 1 asma . 


What  makes  a  MMIS  special  is  that  the  outer  walls  of  the 
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chamber  are  covered  by  a  checkerboard  array  of  permanent 
magnets.  Without  the  magnets  the  electrons  are  quickly 
accelerated  into  the  walls  of  the  chamber  and  lost  from  the 
plasma.  This  results  in  a  low  ionization  efficiency  for  the 
electrons.  With  the  magnets  in  place,  a  complex  multipole 
(or  multicusp  -  hence,  the  name)  magnetic  field  exists  at 
the  wall.  This  field  tends  to  reflect  the  energetic 
electrons  back  into  the  plasma  [2:1822].  Thus  the  ionization 
efficiency  of  the  electrons  is  improved.  This  can  be 
verified  by  looking  at  the  electron  containment  time;  it  is 
found  [3]  that  with  the  magnets  in  place  the  containment  time 
is  much  improved. 

It  has  been  proposed  that  negative  hydrogen  ions  are 
produced  by  these  energetic  electrons  in  a  two-step  process 
[  5:223].  First,  the  energetic  electrons  (E  20  eV)  collide 
with  molecular  hydrogen  and  electronically  excite  the 
hydrogen.  The  excited  hydrogen  then  undergoes  an  electronic- 
vibrational  transition,  producing  vibrational ly  excited 
hydrogen.  This  v ibrat i onal ly  excited  hydrogen  then  collides 
with  slow  electrons  (  E  <_  2  eV)  ,  The  electrons  undergo 
d i ssoc i a t i ve  attachment  onto  the  hydrogen,  producing  H*  . 

Thus  we  have  the  two  processes: 


ei  isi  +  Hj  [X1  2' t  (  v  =  0)  ]  -»  Hi  [B>2»  ,  ,C»Wu  ,  .  .  .  ,V  ]  +  e 

-  Hi  [ X1  2*  g  (  v "  )  ]  +  e  +  hi 


e«  i  o  w  +  Hi  (  v"  )  -»  H*  +  H 


The  negative  ions  produced  can  be  lost  through  several 
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processes  which  include  detachment  induced  by  energetic 
electrons  and  mutual  neutralization  with  positive  ions.  In 
order  to  reduce  losses  caused  by  the  energetic  electrons,  a 
linear  cusp  of  permanent  magnets  is  added  to  the  walls  of  the 
chamber,  as  shown  in  Fig.  1.  These  magnets  introduce  a  sheet 
of  magnetic  field  into  the  plasma  separating  the  discharge 
region  into  two  sections.  The  first  is  called  the  production 
region  and  is  adjacent  to  the  filaments.  In  this  region,  the 
energetic  electrons  produce  the  vibrationaily  excited  H2  . 

The  neutral  species  and  the  thermal  electrons  then  drift 
across  the  magnetic  field  into  the  extraction  region.  It  is 
here  that  the  dissociative  attachment  occurs.  The  energetic 
electrons  are  unable  to  pass  through  the  magnetic  field; 
hence  the  naming  of  this  field  as  the  filter  field.  A  MMIS 
with  such  a  filter  field  is  often  referred  to  as  a  tandem 
d i scharge . 

Further  optimization  of  a  MMIS  is  the  subject  of  several 
studies  that  have  been  performed  [2,6-9].  In  particular, 
much  work  has  been  concentrated  on  consistently  modeling  the 
plasma  discharge.  Bretagne,  et  <3l ,  [10]  numerically  solved 

the  Boltzmann  equation  to  calculate  the  electron  energy 
distribution  function  in  the  H2  discharge.  This  allowed  them 
to  determine  scaling  laws  for  H-  production  in  terms  of 
discharge  current,  discharge  voltage,  and  neutral  gas 
pressure.  These  detailed  calculations  produced  good 
agreement  with  experimental  data;  however,  they  are  very 
tedious  to  perform. 
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Problem  and  Scope 

The  present  study  is  concerned  with  a  simplified 
approach  to  modeling  the  magnetic  multicusp  Hi  discharge 
using  Maxwellian  moments  of  the  Boltzmann  equation.  A 
preliminary  study  of  this  approach  [11]  has  yielded  results 
consistent  with  both  experimental  data  and  detailed 
theoretical  calculations.  This  thesis  takes  a  more  self- 
consistent  approach  towards  modeling  a  MMIS  using  moment 
equations.  This  approach  will  make  it  easier  to  determine 
scaling  laws  for  H*  production  in  terms  of  the  plasma 
parameters,  thus  allowing  the  optimum  operating  regime  of  the 
MMIS  to  be  determined. 


Approach 

This  research  began  with  a  compilation  of  the  cross- 
section  data  for  the  relevant  processes  within  the  hydrogen 
discharge.  This  data  was  then  transformed  into  average 
collision  rates  (cm3/sec)  assuming  a  Maxwellian  distribution. 
Next,  the  density  and  energy  moments  of  the  Boltzmann 
equation  were  developed  with  consideration  to  the  plasma 
chemistry.  These  t ime -dependent ,  coupled  equations  were  then 
solved  numerically,  yielding  the  time  evolution  of  the 
electron  temperature  and  the  densities  of  the  electrons  and 
ions.  The  resulting  data  were  compared  to  experimental  and 
theoretical  data  for  verification  of  the  model.  Finally, 
scaling  laws  with  discharge  parameters  were  determined. 
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CHAPTER  II  -  THE  MOMENT  APPROACH 


In  this  chapter,  the  reactive  and  non-reactive  kinetics 
of  the  hydrogen  discharge  will  be  presented.  Based  on  the 
identified  processes,  Maxwellian  moments  of  the  Boltzmann 
equation  will  then  be  developed.  Finally,  the  numerical 
solution  of  the  moment  equations  will  be  examined. 

Processes  Occurring  in  a  Hydrogen  Discharge 

In  considering  the  relevant  reactions  occurring  in  the 
hydrogen  discharge,  it  is  necessary  to  look  at  production  and 
loss  mechanisms  for  all  of  the  species  being  considered.  In 
this  study  these  species  were:  electrons  (e),  negative 
hydrogen  ions  ( H~ ) ,  positive  hydrogen  ions  (both  H2 ♦  and  H*  ) , 
and  vibrationally  excited  hydrogen  molecules  (Hj(v)).  In 
addition,  the  electron  temperature  was  modeled. 

All  of  the  reactions  will  be  characterized  by  a  rate 
coefficient,  Rj  (E)  (in  cm3 /sec),  that  indicates  the  rate  at 
which  the  reaction  proceeds  for  an  interaction  energy  E. 

Reaction  rates  are  given  by  the  expression 

Ri  (  E )  =  a,  ( E  )  v  (  E  )  (  3  ) 

where  Ri  is  the  reaction  rate  for  process  i  ,  Oi  is  the  cross- 
section  for  an  interaction  energy  E  for  process  i,  and  v  is 
the  relative  velocity  of  interaction  written  in  terms  of 
energy.  Since  the  reactions  dealt  with  here  are  between 
electrons  and  heavier  particles,  it  is  assumed  that  the 
interaction  energy  is  simply  the  electron  energy. 
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In  using  Maxwellian  moments  of  the  Boltzmann  equation, 
it  will  be  necessary  to  use  an  average  value  for  the  reaction 
rates,  <Ri  >.  This  averaging  is  done  over  a  Maxwell-Boltzmann 
distribution : 

<Ri  >  =  S  Oj  (E)  v  ( E )  fm  (E)  dE  (4) 

n« 

with 


fnB(E)  =  2//TT  n,  ( T*  )-  3  / 2  fE  e-«/T* 


5) 


■  * 


v.> 


where  T*  is  the  electron  temperature  (in  eV)  and  ne  is  the 
electron  density  (in  cm-*). 

The  cross-section  data  used  for  most  of  the  reactions 
considered  in  this  paper  are  experimental.  This  prevented 
the  need  to  make  complex  theoretical  calculations  for  this 
moment  approach.  References  to  the  source  of  cross-section 
data  is  given  in  the  form  [data-XX : XXX ]  after  the  reaction  is 
introduced.  For  those  cases  that  pre-calculated  reaction 
rate  data  is  used,  the  reference  is  in  the  form  [rate- 
XX: XXX]  . 


Production  and  loss  of  H~ 


It  has  been  indicated  [6:1247]  that  the  primary  produc¬ 
tion  mechanism  for  H*  in  a  MMIS  is  through  dissociative 
attachment  of  low  energy  electrons  onto  vibrationally  excited 
Ha.  Another  contributing  process  [12:L22]  is  the  dissocia¬ 
tive  recombination  of  an  electron  onto  Ht +  (in  this  study 
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this  will  be  called  polarized  dissociation  to  differentiate 
it  from  dissociative  recombination  resulting  in  neutral 
products).  These  two  reactions  can  be  written  as: 


Dissociative  Attachment  (att)  [ rate- 1 3 : 9 18 ] 
e  +  Hj  ( v )  ->  H  +  H* 


Polarized  Dissociation  (pdiss)  [data- 1 4 : 1 57  2  ] 
e  +  Hj*  -i  H*  +  H-  l't 


Four  loss  mechanisms  of  H"  are  considered.  The  first 
three  are  collisional  detachment  processes  [15:73].  These 
are  detachment  by  an  electron,  associative  detachment  with 
atomic  hydrogen,  and  detachment  by  atomic  hydrogen.  The 
fourth  process  is  mutual  neutralization  with  H*  .  These 
reactions  can  be  written  as: 


Electron  Detachment  (edet)  [data- 16 : 1 347 ] 


e  +  H-  -»  H  +  2e 


Associative  Detachment  (adet)  [ rate- 1 7 : 23  1  ] 


H  +  H-  -»  e  +  Hi 


Hydrogen  Detachment  (hdet)  [ rate- 1 7 : 232 ] 


H  +  H-  ^  2H  +  e 


Mutual  Neutralization  (mneut)  [data- 18 : 437 , 19 : L370 ] 
H*  +  H-  -»  H  +  H  (11) 


These  six  reactions  constitute  the  production  and  loss 
mechanisms  for  H-  ions  that  will  be  considered  in  this  paper 
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The  cross-sections  for  Eqns .  7,  8,  and  11  are  shown  in  Fig. 

2. 

Production  and  loss  of  electrons 

Looking  at  the  electron  energy  distribution  function 
(EEDF)  for  a  typical  MMIS,  shown  in  Fig.  3  (from  Bretagne,  et 
<al  )  ,  it  is  seen  that  there  are  two  distinct  regions  to  the 
distribution.  For  low  energies,  the  distribution  is  charac¬ 
teristic  of  a  Maxwellian  distribution  of  electrons.  For 
higher  energies,  the  distribution  becomes  non-Maxwell ian . 

This  indicates  that  electrons  in  the  discharge  can  be  divided 
into  two  species:  fast  electrons  (E  >  10  eV)  and  slow 
electrons  (E  <  10  eV).  In  this  study,  the  slow  electrons 
(which  are  dominant)  are  the  species  primarily  considered. 

The  fast  electrons  are  treated  as  a  production  mechanism  for 
the  slow  electrons,  for  vibrational ly  excited  hydrogen,  and 
for  positive  hydrogen  ions.  Handling  the  electrons  in  this 
manner  makes  it  possible  to  use  Maxwellian  moments  of  the 
Boltzmann  equation  and  get  good  results,  even  with  the 
distribution  being  non-Maxwe 1 1 ian .  The  fast  electrons  will 
be  treated  in  a  later  section.  This  section  will  deal  with 
the  slow  electrons  only. 

Slow  electrons  are  produced  through  five  reactions.  The 
first  three  are  the  detachment  processes  listed  in  the 
previous  section  (Eqns.  8-10).  The  other  two  are  ionization 
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reactions : 


Hi  Ionization  (h2ion)  (data-20: 1471 ] 

Hi  +  e  -*  Hi*  +  2e  (12) 

H  Ionization  (hion)  [data-21 : 1948] 

H  ♦  e  -*  H4  ♦  2e  (  13 ) 

In  addition  to  these  reactions,  slow  electrons  are 
produced  when  the  fast  electrons  undergo  ionizing  collisions 
This  production  term  will  be  considered  when  the  fast 
electrons  are  studied. 

Electrons  are  lost  through  five  processes  also.  These 
are  dissociative  attachment  (Eqn.  6),  polarized  dissociation 
(Eqn.  7),  three-body  recombination  of  both  H4  and  Hi*,  and 
dissociative  recombination.  These  last  three  reactions  are 
written: 


Ha  4  Recombination  (h2rec) 

Ha  ♦  +  2e  -»  Ha  +  e  (  14  ) 

H4  Recombination  (hrec) 

H4  +  2e  -*  H  +  e  (  15 ) 

Dissociative  Recombination  (dree)  [data-22 : 24  1  ] 

Ha  4  +  e  H  +  H  (  16  ) 

For  the  three-body  recombination  rates,  the  equation 
given  by  Golant,  et  al ,  [23:68]  was  used.  This  equation  is 

<Rr«c>  =  a  ell  n,  *  8.75  x  10**T  T.  • 4  ■  *  n.  (17) 
/nw  T«  4  •  5 
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where  d  is  a  coefficient  of  the  order  of  unity,  e  is  the 


electron  charge,  and  m«  is  the  electron  mass. 

The  cross-sections  for  'qns  .  12,  13,  and  16  are  shown  in 
Fig.  4. 


Production  and  loss  of  Hi  * 


Positive  hydrogen  ions  are  produced  through  two  ioniza¬ 
tion  processes.  The  first  is  through  ionization  (Eqn.  12)  by 
the  slow  electrons.  The  second  is  through  ionization  caused 
by  the  fast  electrons.  This  second  process  will  be  consid¬ 
ered  in  the  section  where  fast  electrons  are  dealt  with. 

The  loss  processes  considered  for  H* *  are  dissociative 
recombination  (Eqn.  16),  polarized  dissociation  (Eqn.  7),  and 
three-body  recombination  (Eqn.  14). 

Production  and  loss  of  H* 


Positive  ions  of  atomic  hydrogen  are  produced  through 
the  processes  of  ionization  (Eqn.  13)  and  polarized  dissocia¬ 
tion  (Eqn.  7).  They  are  lost  through  three-body  recombina¬ 
tion  (Eqn.  15)  and  mutual  neutralization  (Eqn.  11). 


P roduction  and  loss  of  vibrational ly  excited  Ha 


Vibrationally  excited  hydrogen  molecules  are  produced 
through  E-V  processes,  as  was  indicated  in  Chapter  I.  It  has 
been  determined  through  isotope  studies  [6:1247]  that  vibra¬ 
tional  levels  with  6  v  ^  10  are  the  primary  contributors 
to  the  production  of  H-  through  dissociative  attachment  in  a 
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MMIS.  These  levels  are  populated  by  fast  electrons  electroni¬ 
cally  exciting  Hj  and  then  the  hydrogen  undergoing  radiative 
relaxation  into  a  vibrational ly  excited  state  (Eqn.  1) 
[24:51].  A  typical  vibrational  distribution  is  shown  in  Fig. 
5.  It  is  seen  that  these  levels  lie  along  a  plateau  of  the 
distribution,  indicating  that  the  vibrational  densities  for 
levels  6  to  10  are  approximately  equal.  If  the  average 
reaction  rates  for  dissociative  attachment  are  now  looked  at 
(see  Fig.  6),  it  is  seen  that  the  reaction  rate  increases  for 
increasing  vibrational  level.  Combining  a  knowledge  of  the 
vibrational  distribution  with  the  average  reaction  rates  for 
dissociative  attachment  it  is  apparent  that  only  the  vibra¬ 
tional  levels  between  6  and  10  should  contribute,  as  was 
indicated . 

For  this  work  a  representative  population  H2(v),  which 
represents  an  average  value  for  the  vibrational  densities 
between  v=6  and  v=10,  is  used.  In  addition,  for  the  average 
reaction  rate  for  dissociative  attachment,  the  average  for 
the  reaction  rates  between  v=6  and  v=10  is  used. 

The  electronic  excitations  used  to  populate  this  average 


vibrational  level  are: 

e  +  Ha(X)  -  e  +  H2<B’»J*U)  [data-25 : 151 1 ]  (18) 

e  +  H2(X)  -  e  +  H2(B"»Z*u)  [  da  La-25  :  1  5  1  1  ]  (19) 

e  +  H2(X)  -*  e  +  H2(D‘1TU)  [data-25  :  15  1 1  ]  (20) 

e  +  H2(X)  -»  e  +  H2(D’*Tr„)  [  data- 25  :  1  5  1  1  )  (21) 

e  +  H2(X)  -  e  +  H2  (E-F»  !♦  «  )  [  data-25  :  1 5 1  1]  (22) 
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Hj  (X) 

-*  e 

♦  H2  (a3  2*  t  ) 

[data-26 : 1886 ] 

(23) 

e 

+ 

Hj  (X) 

-*  e 

+  H2  (b3Z*  .  ) 

[data-26: 1886] 

(24) 

e 

+ 

H*  (X) 

->  e 

+  Hj  (c3  ir« ) 

[data-26 : 1886 ] 

(25) 

It  is  assumed  in  this  work  that  each  of  these  excitations 
will  proceed  to  relax  into  a  vibrational  state,  producing 
H2 ( v  )  .  The  cross-sections  for  these  reactions  are  given  in 
Fig.  7.  The  choice  of  which  electronic  excitations  .0  use 
was  guided  by  the  work  of  Bretagne,  et  al  [10:814]. 

Vibrat ional ly  excited  hydrogen  is  considered  to  be  lost 
through  three  processes  in  this  study.  It  is  lost  through 
dissociative  attachment,  through  wall  collisions,  and  through 
vibrat ional-translational  (V-T)  relaxation. 

V-T  relaxation  is  given  by  the  reaction  [27:4] 

H2  ( v )  +  H  -+  H 2  (v  =  0)  +  H  (26) 

The  rate  coefficient  for  the  v=l  ->  \  ^0  reaction  is  given  by 
[27:4] 


<Rvt>  S  1.5  x  10'10  e"  ■  1  ®  0  3  T  <  cm3 / sec  (27) 

where  T*  is  the  gas  temperature  in  eV.  Corse,  et  al,  say 
that  the  rate  coefficients  for  higher  vibrational  levels  have 
a  v3  dependence  at  low  gas  temperature,  and  are  linear  with  v 
at  high  gas  temperature.  In  this  study,  a  gas  temperature  of 
about  450  K  was  used,  thus  the  linear  dependence  of  <Rvt>  was 
assumed.  Since  levels  between  6  and  10  are  considered,  <Rvr> 
given  in  Eqn.  27  was  multiplied  by  the  average  value  of  v, 
i.e.-  by  a  factor  of  8. 
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For  wall  losses  of  Hi(v),  it  was  assumed  that  10 
collisions  with  the  walls  were  necessary  to  bring  the  excited 
hydrogen  back  to  the  v=0  state.  The  frequency  of  wall 
collisions  can  be  written  (assuming  free-flight  to  the  wall) 
as  the  ratio  of  Hi  velocity  to  the  average  distance  to  the 
wall,  or,  in  terms  of  the  gas  temperature, 

vvibwiii  =  1.56  x  10#  /T  t  I'l  (28) 

L 

where  L  is  the  average  distance  to  the  walls  of  the 
discharge.  Thus  the  frequency  of  loss  of  v ibrat iona 1 1 y 
excited  hydrogen  to  the  walls  is  vvibwni/10. 

The  fast  electrons 

It  has  already  been  indicated  that  the  electron  distri¬ 
bution  can  be  divided  into  two  distinct  portions.  The  low 
energy  part  of  the  distribution  can  be  approximated  by  a 
Maxwellian  distribution  of  electrons.  In  this  part  of  the 
distribution,  the  electrons  are  dominated  by  electron-elec¬ 
tron  collisions  and  a  thermal  distribution  exists.  The  high 
energy  portion  of  the  distribution  is  where  electron-neutral 
collisions  dominate. 

The  problem  is  in  how  to  approximate  the  high  energy 
tail  of  the  distribution.  These  high  energy  electrons  are 
produced  by  electrons  being  thermionical ly  emitted  from  the 
filaments  in  the  discharge.  These  emitted  electrons  then 
undergo  ionizing  collisions  in  the  plasma  producing  high 
energy,  or  secondary,  electrons,  and  positive  hydrogen  ions. 
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The  strength  of  the  source  (number  of  secondary  electrons 
produced  per  second  per  unit  volume)  is  given  by 

So  =  I  _V (29) 

e  V.i  2  Ei 

where  I  is  the  discharge  current,  V0 1  is  the  plasma  volume,  V 
is  the  discharge  voltage,  and  Ei  is  the  ionization  energy  for 
H2 ( in  Volts).  It  is  assumed  that  each  emitted  electron  will 
have  an  ionizing  collision  with  H2  ,  producing  H2 *  and  a 
secondary  electron.  With  this  assumption,  the  classical 
Thomson  theory  of  ionization  is  used  to  determine  the 
distribution  of  secondary  electrons.  According  to  this 
theory  (see  Appendix  A)  the  secondary  electron  distribution 
is  given  by 


ft  b  o  ■  ( E )  =  Ei  (  1  +  Ei  /V)  0  <  E  <  V  (30) 

<E  +  Ei  )* 

To  determine  the  actual  tail  of  the  EEDF,  a  balance  between 
the  source  production  and  collisional  and  wall  losses  is 
considered.  If  ft.ii (E)  is  the  distribution  of  fast  elec¬ 
trons,  i/c  o  1  1  (E)  is  the  rate  of  electron  impact  collision,  and 
vwai  1  (E)  is  the  rate  of  loss  of  electrons  to  the  walls,  then 
we  have  the  balance 

So  ft  ho.(E)  +  IT’  (E)  =  (31) 

Vcoll(E)  f  t  a  i  I  (  E  )  +  Vw  all  (E)  f  t  a  i  1  (E) 
or 

ft.il  (E)  =  So  ft  h  0  .  ( E )  +  IT’  (E)  0  <  E  <  V-Ei  (32) 

lc  0  I  I  (E)  +  Vw  .  |  |  (  E  ) 
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Vcoi  i  (E)  =  no  2|  [ffi  (E)  v(E)  ]  (33) 

1T*(E)  =  n,  Zi  [ffi  (E+Ei  )  v(E+Ej  )  f  t  .  t  i  (B+Ei  )]  (34) 


and  [10:813] 

vw.n(E)  =  v ( E )  A/Vo  i  (1  -  VP/E)/4  (35) 

where  n0  is  the  neutral  number  density,  the  summation  is  over 
all  of  the  electron  collisional  processes,  v(E)  is  the 
electron  velocity  in  terms  of  energy,  Ei  is  the  threshold 
energy  for  process  i,  IT’  represents  the  production  of  a  fast 
electron  at  an  energy  of  E  due  to  losses  from  the  distribu¬ 
tion  at  higher  energies,  A  is  the  effective  wall  Iosb  area  of 
the  MMIS  chamber  (which  takes  into  account  the  magnetic 
shielding  of  the  walls),  and  Vp  is  the  plasma  potential. 

This  is  the  distribution  that  will  be  used  in  this  study  for 
the  tail  of  the  EEDF. 

If  the  probability  of  an  electron  not  undergoing  a  wall 
collision  before  it  has  the  chance  to  undergo  any  other  type 
of  collision  is  given  by  [11:22] 

P(  E  )  =  1  -  vw.i  i  (E) _  (36) 

Vc  o  I  1  (  E  )  +  Vw  all  (  E  ) 


Maxwellian  distribution  and  the  modified  Thomson  distribu¬ 
tion.  This  superposition  comes  from  being  able  to  consider 
the  electrons  as  composing  two  distinct  species.  In  addi¬ 
tion,  since  number  density  distributions  and  not  probability 
distributions  are  being  used,  the  algebraic  addition  of  the 
two  distributions  will  give  us  the  total  distribution: 

f  (E)  =  fNa  (E)  +  ft  .m  (E)  (38) 

This  distribution  is  graphed  in  Fig.  8  (same  conditions  as  in 
Fig  3).  On  comparison  to  Fig.  3,  it  is  seen  that  approxi¬ 
mating  the  tail  of  the  distribution  as  a  modified  Thomson 
distribution  does  a  good  job  of  representing  the  EEDF . 

For  most  of  the  calculations  (e.g.,  the  rate  coeffi¬ 
cients),  only  the  Maxwellian  portion  of  the  distribution  will 
be  considered.  As  can  be  seen  from  Fig.  8,  the  tail  of  the 
distribution  is  several  orders  of  magnitude  lower  than  the 
Maxwellian  portion.  The  tail  of  the  distribution  becomes 
important  for  the  E-V  processes  that  produce  Ha(v),  since  it 
is  the  high  energy  electrons  that  need  to  be  considered. 

Thus,  for  the  electronic  excitations  producing  Hi(v)  we  have 

<R*  i  *  c  >  =  F  Oj  (E)  v  ( E  )  ft  a.  i  (E)  dE  (39) 

J  ft  .i  i  (E)  dE 

We  have  now  dealt  with  the  production  and  loss  of  all  of 
the  species  that  are  relevant  to  this  study.  In  addition,  we 
have  considered  how  the  fast  electrons  will  enter  into  the 
calculations.  At  this  point  we  will  discuss  the  moment 
equations  that  describe  a  MMIS. 
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The  Maxwellian  Momenta  of  the  Boltzmann  Equation 


In  this  study  we  are  concerned  with  modeling  the 
electron  density,  the  ionic  densities,  and  the  electron 
temperature.  To  do  this,  the  density  moments  for  the 
electrons  and  ions  and  the  energy  moment  for  the  electrons 
will  be  needed.  First  we  start  with  the  basic  moment 
equations  of  the  Boltzmann  equation: 


Density  Moment  (Continuity  Equation)  [31:28] 

8N  +  fl_(N  vi  )  =  Sc  i  ( 

at  8xi 


Energy  Moment  [31:42] 

8_(  1/2  N  m  Vj  Vj  )  +  3_(  1/2  N  m  <Wj  wj  >  )  = 

at  at 

-  8__(  1/2  N  m  v-,  <Wj  Wj  >  )  -  8__(l/2  N  m  Vi  Vj  Vj  ) 

dxj  dxi 

-  3J  1/2  N  m  <wi  wj  wj  >  )  -  a_(N  M  <Wj  wj  >  vj  ) 

ax,  ax, 

+  N  <Fi  Ui  >  +  Ecoi  i  (41) 


where  N  is  the  species  density,  vi  is  the  average  directed 
velocity  of  the  species,  m  is  the  mass  of  the  species,  wi  is 
the  thermal  velocity  of  the  species,  ui  =  Vi  +  Wj  ,  Fi  is 
the  external  force  acting  on  the  species,  Scoi  i  and  Ec o  i  i 
represent  the  change  in  density  and  energy,  respectively,  due 
to  collisions,  angled  brackets  represent  an  average  over  the 
distribution,  and  summation  over  repeated  indices  is  assumed. 
For  this  work  a  Maxwellian  distribution  is  assumed.  The 
densities  are  assumed  to  be  uniform  and  isotropic  over  the 
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plasma  volume.  For  a  MMIS  this  uniformity  is  achieved  by  the 
magnets  shielding  the  walls.  It  is  also  assumed  that  Fi  =  0. 
We  are  considering  the  source  electrons  to  be  part  of  Sc o i i 
and  Ec o i  i  instead  of  producing  an  electric  field.  Also,  the 
magnetic  field  on  the  walls  is  considered  to  only  effect  the 
area  of  possible  wall  loss  of  electrons. 

Consider  first  the  continuity  equation.  For  an  iso¬ 
tropic  distribution  vi  =  0.  Therefore  the  continuity 
equation  becomes 

0N  =  Scoi  i  (42) 

at 

Sc « i  i  for  the  various  species  is  all  that  needs  to  be 
determined  for  the  density  moments.  The  continuity  equations 
for  all  of  the  species  in  the  plasma  can  now  be  written: 


an- 

at 

=  Sc  0  II 

,  I  - 

(43) 

8ne 

at 

=  Scull 

,  • 

(44) 

3n« 

at 

-  Sc  0  1  1 

.Hi 

(45) 

3ne  * 
at 

=  Scot 

1  ,  I  ♦ 

(46  ) 

3nv 

at 

=  Scot! 

.  1  2  (  »  ) 

(47) 

where  n-  is  the  H~  density,  n»  is  the  Hj  *  density,  m»  is  the 
H*  density,  and  n»  is  the  Hj(v)  density. 

Since  the  distribution  of  slow  electrons  is  Maxwellian, 
the  energy  equation  can  also  be  much  simplified.  Note  that 


1/2  m  <wi  wi  >  =  3/2  T« .  With  this  substitution  Eqn.  41  can 
be  written  {still  noting  that  vi  =  0) 

8_(N  3/2  T.  )  =  -  a_(  1/2  N  m  <wi  w,  w;  >  )  Ec  •  n  (48) 

at  at 

To  simplify  the  equation  further  the  first  term  on  the  right 
hand  side  needs  to  be  looked  at.  This  term  is 

<Wj  wj  wj  >  =  1/N  J1  wi  wj  wj  f*»(ui  )  dui  (49) 

f m b ( u i  )  dui  is  an  even  function  of  Wi  ,  therefore  the 
integrand  of  Eqn.  49  is  seen  to  be  an  odd  function  of  Wi  and 
Eqn.  49  is  identically  equal  to  zero  for  a  Maxwellian 
distribution.  Finally,  we  have  (letting  <E>  =  3/2  T» ) 

a_(  n«  <E>)  =  Ecu  (50) 

at 

As  in  the  density  equation,  only  Ec » i  i  needs  to  be  determined 
in  order  to  solve  the  energy  equation. 

The  densi ty  equation  collision  terms 

Looking  at  the  plasma  reactions  involving  H~  we  see  that 
the  collision  term  for  negative  ions  is  given  by  (the  angle 
brackets  on  the  rate  coefficients  are  left  off  for  ease  of 
express i on  ) 

Sc  o  i  i  ,  ■  -  -  n#  n  v  R»  t  i  +  n«  n<  Rp  d  i  ■  ■  —  na,n-R«n*at 

—  n,  n<  R,  ^ ,  t  —  mn.Ridct  ~  nin-Rhd«t  (51) 

For  electrons  we  have 
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Sc  o  i  i  ,  i  -  So  (  ♦  n*  no  Rk  1 1  •  •  +  n*  ni  Rh  i  •  ■  ~  n*  n»  R*  t  t 

-  no n« Rp  dii*  -  n*  n*  Rh  2  r • c  ~  no  n§ ♦  Rh  roc 

+  no  n- Ro  dot  +  n-  ni Ra  dot  —  no  n«  Rd  roc  +  ni n- Rh  dot 

(52) 

Where  <  P>  is  the  average  probability  of  a  fast  electron  not 
undergoing  a  wall  collision  before  it  has  a  chance  to  undergo 
any  impact  collision.  It  is  defined  by 

<P>  =  /  P(E)  ft  ho.(E)  dE  (53) 

where  P(  E )  was  defined  in  Eqn.  36,  and  fthoa(E)  was  defined 
in  Eqn.  30. 

The  collision  term  for  H2 *  is 

Sc  oi  1  ,  in  -  So  <  P>  +  no  no  Rh  2 t  •  0  -  no  n*  Rd  roc 

-  no  n»  Rp  d  i  a  *  ~  no  n,  Rh  aroc  (54) 

The  collision  term  for  H*  is 

Sc  0  1  1  ,  b  ♦  -  no  ni  Rh  i  o  a  t  no  n,  Rp  di  1 1  —  ns  ♦  no  Ra  roc 

~  ni 1 n>  Rodo  u  t  (55) 

The  collision  term  for  H2(v)  is 

Sc  0  1  1  ,  a  2  (  v )  -  nf  n0  2  Rotoc  —  n,  Hv  R,  t  t  —  n,  v,  j  b  wi  i  1  /  1  0 

-  n«  n, Rv  t  ( 56  ) 

where  Ro  1  c c  is  given  by  Eqn.  39,  and  nr  is  given  by 

nf  =  J  ft  ..  1  (E)  dE  (57) 
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and  represents  the  total  number  of  fast  electrons. 


The  energy  equation  collision  term 

The  collision  term  of  the  energy  equation  must  take  into 
account  everything  that  will  cause  the  slow  electrons  to  gain 
or  lose  energy.  The  only  energy  production  mechanism  in  the 
system  is  the  source.  Superelastic  collisions  were  ignored 
in  this  study  since  several  sources  [32:1203,11:22]  indicated 
that  they  have  little  effect  on  the  average  energy  and 
electron  density.  Energy  is  lost  through  electronic  excita¬ 
tions,  through  vibrational  excitations,  through  ionization, 
in  recombination  processes,  in  detachment  processes,  and 
through  momentum  transfer.  Thus  the  collision  term  in  the 
energy  equation  can  be  written  as 

Ecoii  =  Py  -  n,  no  2  E*i«cRei«c  -  n,  n,  E»  i  bR»i  b 

—  n,  no  Eb  i  i  o  a  Rt  >  i  o  i  ~  n«  ni  Eb  i  o  ■  Hh  i  o  i 

—  n,ni <  E  Rd  r  «  c  ( E ) >  —  n, n*  Rh2r*c(E)i 

—  n*  ng  *  <  E  Rhr«c(E)i  —  ne  n» (E  Ratt(E)) 

—  n,n-  <E  R*d*t(E)>  -  n-ng<E  Rhd*«(E)> 

—  n-  ng  <  E  R,  om  (  E  )  >  -  2  (  m,  /M  )  n,  n0  <  E  R.t(E)> 

(58) 

where  <E9 >  is  the  average  energy  of  the  source  and  is  defined 
by 

<  Eg  >  =  r  E  P(E)  ft  ho  .  (E)  dE  (59) 

J*  P(  E)  ft  ho.(E)  dE 

Ei  is  the  threshold  energy  for  process  i  and  M  is  the  mass  of 
H2  .  Th  e  subscript  vib  refers  to  the  excitation  of  Hj  from 
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the  v=0  state  to  the  v=l  state.  The  subscript  at  refers  to 
the  momentum  transfer  reaction.  For  the  energy  equation, 

R«  1  •  c  is  averaged  over  the  Maxwell-Boltzmann  distribution, 
and  not  the  modified  Thomson  distribution,  since  the  slow 
electrons  are  being  considered. 

It  should  be  noted  that  in  Eqn.  58  the  recombination  and 
detachment  terms  are  averages  of  the  rate  coefficients 
multiplied  by  the  electron  energy.  As  a  simplification  for 
this  model,  this  average  will  be  approximated  by  the  average 
electron  energy  times  the  average  reaction  rate.  This  will 
result  in  a  slightly  less  accurate  electron  energy,  but  it 
will  reduce  the  computations  needed.  Thus  the  energy  equa¬ 
tion  collision  term  actually  used  in  the  calculations  in  this 
study  is 

Ecol  I  -  <  E»  >  So  ^  P)  -  ne  no  2  EelocRelec  —  ne  no  Ev  i  b  Rv  i  b 

-  n,  no  Eh2  i  o  iRbi  i  o  b  -  n,  niEti  o  iRbi  o  i 

-  n,  <EHRd  t  »  c  (  E)  )  -  n*  n«  (E>  (Rb  j  r  t  c  (E)i 

-  n«nH*  <EXRbr.c  (E)>  -  n»  n»  <E>  <Ri  t  t  (  E)  > 

-  n*  n- <E>  <R«  d <  t  ( E ) >  -  n- nH <E> <Rn d » t  (  E )  > 

-  n-  nH  <E>  <R»  d  e  t  (  E  )  >  —  2(mo/M)n*no^E^R»t  (60) 

Solution  of  the  Coupled  Moment  Equations 

An  analytic  solution  of  the  moment  equations  is  not 
possible  due  to  the  experimental  nature  of  the  cross-section 
data.  It  was,  therefore,  necessary  to  solve  the  equations 


numerically.  A  complete  description  of  the  numerical  tech- 


niques  used  and  a  listing  of  the  actual  computer  code  are 
included  in  Appendix  B.  In  this  section  only  a  brief  summary 
of  the  solution  process  will  be  presented. 

The  equations  are  solved  by  using  a  Runge-Kutta-Fehlberg 
method  of  orders  4  and  5  [29:129-147].  This  is  an  adaptive 
routine  and  is  stable  for  mildly  stiff  systems  of  differen¬ 
tial  equations.  The  average  rate  coefficients  as  a  function 
of  the  average  electron  energy  are  stored  in  the  program  in 
look-up  tables.  These  tables  were  calculated  using  cross- 
section  data  and  integrating  (using  an  adaptive  quadrature 
routine  [30:174-175])  over  either  f*»  or  ft»ii  . 

The  program  proceeds  until  a  steady  state 
(mathematically  determined  by  the  time  derivatives  in  the 
moment  equations  equaling  zero)  is  determined  to  be  achieved 
or  until  a  set  time  is  reached  (typically  20  msec),  whichever 
comes  first. 

In  the  next  chapter  the  results  of  this  model  are 
compared  to  various  experimental  and  theoretical  calcula- 


t  i  ons  . 


CHAPTER  III  -  RESULTS  AND  ANALYSIS 


In  this  chapter  the  results  of  the  model  presented  in 
the  previous  chapter  will  be  discussed  and  compared  to  other 
experimental  and  theoretical  work.  Scaling  laws  with  the 
discharge  parameters  will  then  be  presented.  These  laws  will 
then  be  used  to  estimate  the  best  operating  conditions  for  a 
MMIS  . 

Comparison  of  Model  to  Previous  Work 

For  a  preliminary  consistency  check  of  the  moment  model, 
the  output  was  compared  to  the  theoretical  calculations  of 
Bretagne,  et  alt  [10],  and  the  experimental  data  of  P6alat, 
et  alt  [33]  for  a  pressure  of  40  mTorr,  a  discharge  voltage 
of  90  V,  a  plasma  potential  of  2  V,  and  56%  dissociation  of 
H2  into  H.  The  discharge  currents  considered  were  1-10  A. 
Table  I  presents  the  results  obtained. 

The  electron  temperature  and  electron  density 
calculations  are  very  good  for  the  current  range  considered. 
In  fact,  they  are  generally  better  predictions  of  the 
experimental  data  than  Bretagne’s  work.  The  ratio  of  n-/n», 
on  the  other  hand,  is  as  much  as  an  order  of  magnitude  off 
from  the  experimental  data.  This  is  probably  due  to  my 
treatment  of  the  vibrational  population  density.  The  use  of 
a  representative  population,  instead  of  using  the  full  set  of 
vibrational  master  equations,  causes  a  loss  of  information. 
This  loss  of  information  results  in  the  negative  ion  density 
calculation  being  an  order  of  magnitude  estimate. 


TABLE  I 


Comparison  of  Discharge  Parameters 

40  mTorr,  90  V,  56%  dissociation,  Vp  =  2  V 
T <  =  450  K,  Vo  i  =  10  1,  A  =  700  cm*,  L  =  14.3  cm 

(a)  Present  work  (b)  Experiment  [33]  (c)  Theory  [10] 


1(A) 

T, (eV) 

n«  ( 101 1  cm- 

3  ) 

n-  /nt  ( 10- 

3  ) 

1 

.  43* 

.  43b 

.  3 2C 

2.5* 

1  .  7b 

1 . 9C 

.65*  7.0b 

7 . 0C 

3 

.50 

.65 

.42 

4.1 

5.7 

3 . 8 

1.8  4.0 

4  .  1 

10 

.62 

.85 

.  59 

7.1 

9.8 

6 . 7 

4.8  - 

2.7 

TABLE  II 


Comparison  at 

Low  Pressure 

2  mTorr,  50 
T *  =  450  K, 

V,  5  A,  .4% 
V0i  =  8.8  1, 

dissociation,  VP 
A  =  830  cm2  ,  L 

=  1.15  V 
=  1  4 . 3  cm 

(  a ) 

Present  work 

(b)  Experiment  [24]  (c) 

Theory  [24] 

n*  (  1  0 1  0  cm-  3 

)  T.  (eV) 

n-  (10®  cm- 3  ) 

n-  /n«  (  10- 3  ) 

( a  ) 

26 . 30 

.748 

20. 12 

76 . 52 

(b) 

1  .  1-3  .  4 

.75 

3 . 2 

94 . 1-290.9 

( c  ) 

3  .  43 

.  40 

2 . 3 

67.1 
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It  is  seen  from  Table  I  that  the  present  model  predicts 
that  the  ratio  of  n-  /n#  will  increase  with  current,  whereas 
both  Pealat  and  Bretagne  show  a  decrease  with  current.  This 
discrepancy  can  be  partially  accounted  for  by  the  fact  that  a 
dissociation  of  56%  was  used  at  all  currents.  It  will  be 
shown  later  that  the  negative  ion  density  is  inversely 
proportional  to  the  atomic  hydrogen  density.  As  the  current 
is  increased,  you  have  a  larger  electron  density  with  which 
H2  can  be  dissociated.  Therefore,  at  higher  currents,  the 
dissociation  should  be  larger  and  the  total  negative  ion 
density  should  decrease.  For  a  fixed  dissociation,  however, 
it  will  later  be  seen  that  the  behavior  this  model  predicts 
has  been  experimentally  observed. 

The  model  was  also  compared  to  the  work  of  Gorse,  et  al, 
[24]  at  a  pressure  of  2  mTorr,  a  voltage  of  50  V,  a  current 
of  5  A,  a  plasma  potential  of  1.15  V,  and  a  dissociation  of 
.4%.  The  results  of  this  comparison  are  presented  in  Table 
II  . 

The  moment  model  has  predicted  the  electron  temperature 
extremely  well.  The  electron  density  and  negative  ion 
density  calculations,  however,  are  both  off  by  an  order  of 
magnitude.  The  error  in  the  electron  density  calculation  is 
most  likely  due  to  the  fact  that  wall-losses  were  ignored  in 
the  electron  continuity  equation.  For  the  data  in  Table  II, 
a  plasma  potential  of  1.15  V  was  used,  compared  to  2  V  for 
the  data  in  Table  I.  At  this  lower  plasma  potential,  more  of 
the  slow  electrons  can  overcome  the  potential  barrier  and  be 


lost  to  the  walls.  This  increased  wall-loss  could  possibly 
explain  why  the  observed  electron  density  is  lower  than 
calculated  by  this  model.  The  negative  ion  density  being  an 
order  of  magnitude  too  high  is  then  due  directly  to  n#  being 
an  order  of  magnitude  too  high,  since  n-  0[  n«  . 

With  fairly  good  success  of  the  model  for  various  values 
of  the  discharge  parameters,  the  code  was  then  run  over  the 
pressure,  voltage,  and  current  parameter  spaces.  The  results 
of  these  runs  were  used  to  determine  scaling  laws. 

Scaling  Laws  with  Discharge  Parameters 

The  scalings  of  the  various  particle  species  are 
determined  by  looking  at  the  steady  state  solutions  of  the 
moment  equations.  These  solutions  are  given  by  the  following 
coupled  set  of  equations: 

n-  = _ n,n,Rlt  t  +  n,  n,  RP  a  .  »  » _  (60) 

ni,Rin«iit  +  n,R,d,  t  +  ni  (  Ri  d  •  t  ^  Rh  a  *  t  ) 


n,  =  S9<P>  +  n-  ng  R«  4 «  t  ♦  n-ngRh<let 

X 

(61) 

X  —  nvRatt  +  rw  (  Rp  <J  i  •  s  +  Rh2rec  +  Rdrec) 

~  rio  Rb  2  i  ob  “  Hg  Rb  i  on  ~  11-  R#  d  e  t 

+  ng  ♦  Rb  r  »  c 

m  =  So  <P>  +  no  Rb  2  i  o  a 

n«(Rdr*c  +  Rb  2  r  *  c  +  Rpdisi) 

(62) 

na  ,  =  n«  (  na  Rb  i  on  +  n  *  Rd  d  i  a  *  ) 

n«Rbrec  +  H-Rtne  at 

(63) 

nv  -  Ilf  no  2  Re  i  e  c 

n(  Rj  t  t  +•  naRir  +  vvibwaii/10 

(64) 

S: 


By  looking  at  the  relative  magnitudes  of  the  various 
terms  in  the  above  equations  it  is  possible  to  determine  the 
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dominant  processes  in  the  steady  state.  Using  the  data  from 
running  the  code  between  20-100  V,  1-100  A,  1-100  mTorr,  and 
56%  dissociation,  it  was  found  that  the  densities  are  given 


n-  -  n«  nvRm 

niRittt 


P  >  4  bT 


n*  ~  n* 


/_ S^.<P>  \  1 

ln+  Rd  ric/ 


ng  *  =  n«  ru  Rp  a  t  « » 
n-  Ra seat 

n v  —  _ n t  np  2  Re  i  e  c _ 

neRitt  +  ngRvt  +  Vv  i  b  w  »  i  l  /  1 0 


As  indicated,  Eqn.  65  is  only  valid  if  the  pressure,  P, 
is  greater  than,  or  about,  4  mTorr.  For  pressures  lower  than 
this,  Eqn.  65  will  produce  an  answer  with  about  12%  error. 

To  be  more  accurate,  the  polarized  dissociation  and  mutual 
neutralization  terms  of  Eqn.  60  must  be  retained.  This 
indicates  that  at  lower  pressures,  polarized  dissociation  and 
mutual  neutralization  start  becoming  the  dominant  production 
and  loss  mechanisms  for  negative  hydrogen  ions.  Thus  we  have 


P  <  4  uT 


n i  *  Rim 


If  the  current  is  greater  than  40  A  and  the  pressure  is 
less  than  10  mTorr,  then  Eqn.  68  can  be  simplified  to 


_ _ nf  np  I  R«  i  n _ _ 

fl(  Ru  t  +  Vvibwall/10 


I  >  40  A  (70) 

P  <  10  mT 


Using  Eqns .  65-70  it  will  be  possible  to  elucidate 
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We  will  first  determine  how  the  electron  density  scales 
with  current,  using  Eqn.  66.  By  refering  to  Eqns .  53,  36, 

35,  33,  and  30  it  is  seen  that  <  P>  is  independent  of  current. 
Also,  the  rate  of  dissociative  recombination  is  independent 
of  current  if  you  assume  that  the  electron  temperature  is 
constant.  Eqn.  29  shows  that  the  source  strength  is  directly 
proportional  to  the  current;  therefore,  we  have 

n,  fl[  I»/*  (71) 

In  Fig.  9,  a  graph  of  n«  as  a  function  of  current  is 
given.  From  the  graph  it  is  obvious  that  the  scaling  given 
by  Eqn.  71  is  valid  for  the  model.  This  scaling  was  also 
found  by  Bailey  and  Jones  [11:21]. 

To  determine  the  scaling  of  the  vibrational  population 
with  current  we  must  first  see  how  nr  scales  with  current. 

The  fast  electrons  will  scale  the  same  as  ftaii  with  respect 
to  the  discharge  parameters  (see  Eqn.  57).  Refering  to  Eqn. 
32  and  ignoring  the  IT’  term,  it  is  seen  that  ftaii  is 
directly  proportional  to  current.  Likewise,  m  will  be 
directly  proportional  to  the  current.  R« i « c  and  R» t »  are 
both  independent  of  current;  therefore,  using  Eqn.  70  and 
ignoring  the  wall-loss  term,  we  find 

n,  (I  I1'1  (72) 
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Fig.  10  is  a  graph  of  n»  as  a  function  of  current  and  shows 
the  dependence  given  by  Eqn.  72. 

With  a  knowledge  of  how  n.  and  n»  scale  with  current  it 
is  possible  to  determine  how  the  negative  ion  density  scales 
with  current.  Again  assuming  that  rates  are  constant  and 
assuming  that  n§ ♦  is  independent  of  current  (we  will  show 
this  shortly),  both  Eqn.  65  and  Eqn.  69  give 

n-  CE  I  (73) 

This  scaling  is  seen  to  hold  in  Fig.  11,  and  is  verified 
experimentally  by  York,  et  al  [35:682]  and  Bacal,  et  al 
[36:23]  . 

Using  the  scalings  for  the  electron  density  and  the 
negative  ion  density,  Eqn.  67  indicates  that  n« ♦  should  be 
constant  with  respect  to  current  changes.  Fig.  12  is  a 
graph  of  n««  versus  current;  the  density  is  fairly  constant 
for  currents  above  10  A.  Below  10  A,  the  inconsistency  is 
due  to  the  fact  that  the  n« ♦  density  had  not  reached  a  true 
steady  state  when  the  code  output  the  density.  This  will  be 
seen  when  the  time  evolution  of  the  densities  is  considered 
later.  The  rest  of  the  calculations  will  not  be  affected  by 
this  lack  of  steady  state  since  none  of  the  other  densities 
depend  on  n§ .  . 

Using  Eqns .  71  &  73  it  is  possible  to  find  the  scaling 
for  the  ratio  of  negative  ions  to  electrons.  This  ratio  is 
important  in  the  extraction  of  negative  ions  since  a  large 
population  of  electrons  extracted  with  the  negative  ions 


’’i  a  t] 


degrades  the  quality  of  the  ion  source  output  [34:56].  The 
ratio  scales  as 


n-_  a  1 1  /  *  (  74  ) 

n. 

A  graph  of  this  ratio  is  seen  in  Fig.  13,  confirming  the 
given  current  scaling. 

Next,  the  scaling  of  the  densities  with  respect  to 
pressure  needs  to  be  examined. 

Scaling  with  pressure 

If  it  is  assumed  that  <  P>  and  the  rate  coefficients  are 
independent  of  pressure,  then  Eqn .  66  indicates  that  the 
electron  density  should  be  independent  of  pressure.  Fig.  14 
which  is  a  graph  of  n«  as  a  function  of  pressure,  shows  that 
the  electron  density  is  basically  independent  of  pressure. 
The  slight  increase  with  pressure  that  is  seen  is  due  to  the 
fact  that  <P>  is  weakly  dependent  on  pressure  (see  Fig.  15), 
since  vc o  i  i  is  a  function  of  pressure. 

In  order  to  determine  the  pressure  scaling  of  nv ,  we 
have  to  look  at  the  scaling  of  nr  with  pressure.  Eqn.  32 
indicates  that  the  tail  of  the  distribution  scales  as 
(P  +  k)*1 ,  where  k  is  a  constant,  and,  thus,  nf  scales  the 
same.  Since  n0  is  proportional  to  pressure,  and  assuming 
that  n*  is  independent  of  pressure,  we  get 


Fig.  13  -  The  ratio  n- /n*  as  a  function  of 
discharge  current  for  the  same  conditions  as  in 
Fig.  9. 


Fig.  14  -  Electron  density  as  a  function  of  gas 
pressure  for  a  90  V,  1  A,  56%  dissociation 
d  i  scharge . 
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Thus  for  low  pressures  the  vibrational  density  should  be 
roughly  proportional  to  pressure  and  it  should  saturate  at 
higher  pressures.  The  vibrational  density  versus  pressure  is 
shown  in  Fig.  16  and  it  is  seen  to  basically  follow  the 
scaling  of  Eqn.  75. 

With  the  pressure  scalings  for  n*  and  n«  it  is  now  easy 
to  determine  how  the  negative  ions  will  scale  with  pressure 
in  this  model.  Since  m  (I  P  we  have,  using  Eqn.  65, 

n-  d  1  (76) 

P  +  k 

Fig.  17  demonstrates  this  scaling  for  the  negative  ions.  For 
high  pressure  the  density  is  inversely  proportional  to  P,  and 
for  low  pressures  signs  of  saturation  are  evident.  Actually, 
the  curve  in  Fig.  17  has  a  maximum  at  about  2  mTorr.  This 
maximum  has  been  experimentally  observed  by  Bacal,  et  al 
[36:22].  The  maximum  appears  to  occur  from  the  change  in  the 
negative  ion  production  and  loss  mechanisms.  This  is 
represented  mathematically  by  the  change  from  Eqn.  65  to  Eqn. 
69  at  these  low  pressures. 

Using  Eqns .  67  &  76  it  appears  that  the  Hf  density 
should  be  given  by 

n8»  d  P  (77) 

Looking  at  Fig.  18,  this  scaling  does  not  seem  to  hold. 
However,  as  was  mentioned  earlier,  the  H*  density  had  not 
reached  a  steady  state  value  at  the  time  of  output  from  the 
code,  so  the  above  scaling  should  not  be  expected  to  hold  for 
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the  data  presented. 

Finally,  the  ratio  of  negative  ions  to  electrons  is 
expected  to  follow  the  same  scaling  as  for  n-  (Eqn.  76). 
This  is  confirmed  in  Fig.  19. 

Next,  we  will  determine  the  scaling  with  discharge 
voltage . 


v 

I 

V 


Scaling  with  voltage 


The  electron  number  density  voltage  scaling  is  due  to 
the  scaling  of  the  source  strength.  From  Eqn.  29,  it  is  seen 
that  So  is  directly  proportional  to  V.  Thus  we  expect  to 
have 


n«  OC  V1  *  (78) 

The  electron  density  as  a  function  of  discharge  voltage  is 
graphed  in  Fig.  20.  The  curve  in  this  figure  follows  the 
scaling  of  Eqn.  78  only  approximately.  Eqn.  78  should  only 
be  used  as  an  indication  that  the  electron  density  is  a 
slowly  increasing  function  of  discharge  voltage. 

If  we  assume  Eqn.  78  is  valid  in  order  to  analyze  the 
voltage  scaling  of  the  rest  of  the  species,  then  we  need  to 
find  the  scaling  of  nr  with  voltage  to  see  how  the 
vibrational  population  scales.  From  Eqn.  32,  we  find  that 
ft  aii  ,  and  hence,  nt  ,  is  directly  proportional  to  the 
voltage.  Combining  this  with  Eqn.  78,  we  get  that  nv  scales 
as 


n»  I  V*  * 
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Fig.  21  shows  nv  versus  voltage.  This  graph  is  roughly 
represented  by  Eqn .  79. 


Using  Eqn.  65  and  incorporating  Eqns .  78  &  79,  we  find 

that 


n-  I  V 


The  actual  scaling  is  shown  in  Fig.  22.  This  graph  basically 
validates  Eqn.  80.  The  scaling  given  by  this  equation  has 
been  experimentally  shown  to  exist  by  Leung,  et  al  [9:365]. 

In  that  study,  it  was  found  that  the  negative  ion  density  was 
directly  proportional  to  voltage  up  to  around  120  V,  at  which 
point  it  saturates. 

The  scaling  of  the  ratio  n- /n*  is  found  from  Eqns.  78  & 
80  to  be 


n-  AC  V»/ 2 


This  ratio  is  graphed  in  Fig.  23.  As  with  the  rest  of  the 
voltage  scalings,  Eqn.  81  only  approximately  duplicates  the 
behavior  that  the  computer  code  indicates  exists.  Part  of 
the  differences  between  these  equations  and  what  is 
calculated  is  that  <  P>  is  a  slowly  varying  function  of  V  (see 
Eqns .  35  &  30 ) . 

From  Eqn.  67,  we  see  that  the  H*  density  should  be 
independent  of  the  discharge  voltage.  Fig.  24  indicates  that 
this  is  not  the  case.  However,  the  fact  that  the  curve  in 
Fig.  24  is  dependent  on  the  voltage  is  due  to  ns*  not  yet 
reaching  a  steady  state,  as  indicated  earlier. 
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Fig.  22  -  Negative  ion  density  as  a  function  of 
discharge  voltage  for  the  same  conditions  as  in 
Fig.  20. 
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Up  to  this  point,  we  have  discussed  the  scalings  of  the 
species  densities  with  respect  to  the  discharge  parameters. 
Next,  we  will  determine  how  the  electron  temperature  varies 
with  the  discharge  parameters. 

Electron  temperature  scaling 

Using  the  energy  equation  (Eqn.  60)  in  a  steady  state 
and  only  considering  the  dominant  terms,  we  find  that  the 
electron  temperature  is  given  by  the  balance 

( Eg  )  So  ^  —  n,noEfi  bRvi  b  t  n#n+(E^Rdr*c  (82) 

+  2  ( m«  /M )  n«  no  ^  E  )  R«  t 


The  relative  strengths  of  the  terms  on  the  right-hand 
side  of  this  equation  vary  depending  upon  the  specific 
current,  pressure,  and  voltage  range  the  source  is  being  run 
in.  In  general,  for  currents  less  than  10  A  and  pressures 
greater  than  10  mTorr,  the  momentum  transfer  term  will 
dominate  over  the  recombination  term.  For  higher  currents 
and  lower  pressures  the  reverse  is  true.  For  all  values  of 
the  discharge  parameters,  the  vibrational  term  is  the 
largest,  representing  anywhere  from  50%  to  90%  of  the  total 
contribution  from  the  right-hand  side.  Since  this  is  true, 
we  will  ignore  both  the  momentum  transfer  term  and  the 
recombination  term  in  order  to  make  it  easier  to  see  how  the 
electron  temperature  scales  with  the  discharge  parameters. 

In  this  case,  Eqn.  82  reduces  to 

^  E  a  ^  So  ^  Py  —  n*noEvibRvib  (83) 
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We  will  now  make  the  assumption  that 

Rv i  b  *  R.e-‘«/T*  (84) 

This  assumption  should  be  valid  since  the  electron 
temperature  is  in  the  neighborhood  of  Evib.  The  constants  Ro 
and  Eo  are  determined  by  fitting  R»i  b  to  Eqn.  84.  Using  Eqn. 
84  in  Eqn.  83,  we  get 

1  =  In  f  (n>nBEvi  bR0  )/(<E.  >So  <P>  )  1  (85) 

T«  Eo 

We  know  how  all  of  the  terms  in  Eqn.  85  scale  with  the 
discharge  parameters,  so  we  find  that 

1  <E  In  P  (86) 

T.  (I  V)»/2 

Figs.  25-27  show  the  electron  temperature  data  as  functions 
of  current,  pressure,  and  voltage,  respectively.  It  is  seen 
that  Eqn.  86  is  a  good  representation  for  what  is  observed. 


Summary  of  Scaling  Laws 

The  scaling  laws  in  terms  of  pressure,  current,  and 
voltage  are  given  by 


n-  OC  I  V 

(87) 

P  +  k 

n,  OC  (I  V)‘  /  * 

(88) 

n»  OC  (I  V)i  /  *  P 

(89) 

P  +  k 

nn  I  P 

(90) 

v  • 
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Fig.  27  -  T» "  1  as  a  function  of  In  (V)_l  1  for  a  40 
mTorr,  1  A,  56%  dissociation  discharge. 
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where  k  is  a  constant  that  depends  on  the  gas  temperature, 
the  plasma  volume,  the  effective  wall  loss  area,  the  plasma 
potential,  and  the  reaction  rates  for  the  electronic 
excitations.  These  proportionalities,  in  addition  to  Eqn. 

86,  give  a  means  to  optimize  the  production  of  negative 
hydrogen  ions  in  a  MMIS,  with  respect  to  the  discharge 
parameters . 

The  Optimal  Operating  Regime  for  a  MMIS 

Using  the  scaling  laws  predicted  by  this  model,  it  is 
easy  to  determine  for  which  values  of  the  discharge 
parameters  the  efficiency  of  a  magnetic  multicusp  ion  source 
can  be  optimized.  To  maximize  the  production  of  negative 
ions,  the  discharge  parameters  must  be  chosen  so  as  to 
maximize  Eqn.  87.  It  is  obvious  that  this  will  be  for  a  high 
current,  high  voltage,  low  pressure  source.  For  voltage  and 
current,  however,  there  are  optimum  values.  As  it  was 
mentioned  earlier,  the  negative  ion  density  actually 
saturates  for  a  voltage  around  120  V.  This  then  represents 
the  voltage  limit  for  efficient  operation.  The  optimum  value 
for  the  pressure  was  found  to  be  between  2  and  4  mTorr. 


In  addition  to  Eqns .  86-100  being  useful  for  determining 
the  optimum  operating  regime  for  a  generic  MMIS,  they  should 
be  useful  to  the  experimentalist  in  the  optimization  of  his 


Before  concluding,  it  is  worth  mentioning  how  the 
densities  and  electron  temperature  evolve  with  time.  For  a 
40  mTorr,  90  V,  1  A  discharge  most  of  the  species  reach  a 
pseudo-steady  state  within  about  2  msec  (see  Figs.  28-33). 

For  the  H*  density  this  is  not  true  (Fig.  33);  it  has  not 
reached  a  steady  state  after  20  msec  (at  which  point  the  code 
was  stopped).  Since  the  H*  density  never  truly  reached  a 
steady  state,  the  scalings  developed  above  are  not  exactly 
seen  from  the  data.  None  of  the  other  parameters  being 
measured  are  dependent  upon  ng  ♦  ,  so  the  fact  that  it  has  not 
reached  steady  state  does  not  affect  the  other  calculations. 
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Fig.  30  -  Temporal  evolution  of  H> 4  density  for  the 
same  conditions  as  in  Fig.  28.  _ 
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Fig.  31  -  Temporal  evolution  of  vibration ally 

excited  hydrogen  density  for  the  same  conditions  a 
in  Fig.  28 . 
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CHAPTBR  IV  -  CONCLUSIONS  AND  RBCOMMBNDATIONS 


Cone  1 ua i ona 

This  study  has  developed  a  method  to  model  a  magnetic 
mulieusp  ion  source  using  self-consistent  Maxwellian  moments 
of  the  Boltzmann  equation.  The  departure  from  a  Maxwellian 
electron  distribution  was  considered  and  incorporated  into 
the  model  by  modifying  the  high  energy  tail  of  the 
distribution  using  the  classical  theory  of  ionization.  The 
distribution  tail  so  derived  is  found  to  scale  as  the  ratio 
of  current  to  pressure,  as  was  previously  seen  from  the 
complete  solution  of  the  Boltzmann  equation  (10:819-820). 

A  numerical  comparison  of  the  results  of  this  model  to 
previous  experimental  and  theoretical  work  show  that  they  are 
in  fairly  good  agreement.  The  most  useful  result  from  this 
work,  however,  is  the  determination  of  how  the  various 
species’  densities  scale  with  the  discharge  parameters. 

The  scaling  laws  obtained  from  this  model  are  supported 
by  several  experimental  and  theoretical  studies.  In 
particular,  the  experimentally  observed  scaling  of  negative 
hydrogen  ions  with  current,  pressure,  and  voltage  have  been 
predicted  by  this  model.  With  these  scaling  laws  it  has  been 
determined  that  for  the  most  efficient  production  of  negative 
hydrogen  ions,  the  discharge  should  be  run  at  high  currents, 
pressures  around  2-4  mTorr,  and  voltages  around  100-120  V. 

The  soiling  laws  allow  negative  ion  source  users  to  determine 
how  t o  vary  the  discharge  parameters  most  efficiently  to 
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current  of  1 4 0  — 180  mA  with  an  omittance  of  about  .07  cm-arad. 
A  typical  MMIS  [5:234,35:684]  produces  a  current  of  about  1 
mA  with  an  emittance  of  around  .03  cm-mrad.  The  factor  of  2 
difference  in  emittances  indicates  that  the  NPB  from  a  MMIS 
will  have  roughly  twice  the  effective  range  as  a  NPB  from  a 
surface  source.  If  it  is  assumed  that  the  current  delivered 
by  the  Penning  source  is  the  minimum  allowed  for  a  weapon 
system  to  be  considered  lethal,  then  the  beam  delivered  by 
the  MMIS  cannot  be  considered  lethal.  However,  it  has  been 
determined  [ 37  :S79]  that  the  total  beam  energy  (and  hence 
beam  current)  required  for  use  in  discrimination  may  be  a 
factor  of  around  100  less  than  required  for  lethality. 
Therefore,  a  MMIS  is  an  excellent  option  for  use  in  a  NPB 
discrimination  system. 

A  final  factor  in  favor  of  volume  sources  over  surface 
sources  is  the  use  of  cesium  in  the  latter.  The  cesium, 
being  highly  reactive,  will  decrease  the  useful  lifetime  of 
the  surface  source.  A  volume  source  does  not  have  this 
prob 1  cm . 

Recom  me  ndatio ns 

To  improve  the  results  obtained,  several  things  can  be 
done  with  this  model.  First,  all  of  the  runs  were  made  with 
a  dissociation  of  56%  (taken  from  the  value  in  Bretagne,  et 
al  (10]).  Adding  a  self-consistent  dissociation  into  the 
model  by  considering  the  appropriate  dissociation  chemistry 
would  probably  impr.  e  the  numerical  results  obtained 


tremendously . 

Second,  the  plasma  potential  was  always  taken  to  be  2  V 
(except  for  the  2  mTorr,  50  V,  5  A  discharge).  Bretagne,  et 
al,  (32:1206-1207)  have  shown  that  a  variation  in  the  plasma 
potential  dramatically  effects  the  Maxwellian  portion  of  the 
KEDF.  The  present  model  could  be  improved  by  adding  a  self- 
cons  intent  calculation  of  the  plasma  potential  based  on  the 
balance  of  charged  particles  falling  on  the  chamber  walls. 

’!  h  i  h  calculation  would  be  complex  due  to  the  nature  of  the 
multieusp  field  at  the  wall. 

Finally,  this  model  ignored  the  presence  of  Hs *  ,  which 
is  actually  the  dominant  ionic  species  found  in  many 
multicusp  discharges  (7:1802).  It  is  indicated  by  Bruneteau 
(15:381)  that  for  a  4  mTorr  discharge,  H*-Hj‘  mutual 
neutralization  is  the  dominant  loss  mechanism  for  negative 
ions,  being  a  factor  of  10  larger  than  the  losses  due  to 
associative  detachment.  In  addition,  she  indicates  that  Hj * 
ions  play  an  important  role  in  the  production  of 
v i bra t i ona 1 1 y  excited  hydrogen.  Thus,  inclusion  of  the 
plasma  chemistry  for  Hj *  should  improve  the  calculations  of 
this  model  significantly. 

Hie  next  logical  step  would  be  to  expand  the  approach  to 
include  <d.r  spatial  variation  of  the  densities  and  electron 
tempo r a t ure .  This  could  possibly  be  accomplished  by  not  only 
using  the  continuity  equation  and  energy  equation,  but  also 
the  Momentum  moment  of  the  Boltzmann  equation.  With  this 
adapt  at  ion  to  the  model,  it.  would  be  easier  to  analyze  what 
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APPBNDIX  A  -  THE  CLASSICAL  THEORY  OF  IONIZATION 
AND  THE  THOMSON  DISTRIBUTION 

In  this  appendix,  the  classical  theory  of  ionization 
will  be  developed.  This  development  will  lead  to  the 
expression  for  the  Thomson  distribution  of  secondary 
electrons  used  in  Chapter  II.  The  development  used  is  that 
given  by  Thomson  and  Thomson  in  Ref.  28,  pages  96-98  and  200. 

We  will  regard  the  ionization  of  a  molecule  by  an 
electron  as  a  collision  between  two  particles  -  the  electron 
and  one  of  the  electrons  in  the  molecule.  Ionization  occurs 
when  the  molecular  electron  receives  energy  greater  than  or 
equal  to  the  ionization  energy  from  the  incident  electron. 
Although  we  will  be  concerned  with  the  collision  between  two 
electrons,  the  following  development  i3  for  any  incident 
charged  particle. 

Let  M  be  the  mass  of  the  incident  particle  and  let  E  be 
its  charge.  Assume  that  its  velocity  before  the  collision  is 
Vo  and  let  b  be  the  impact  parameter  for  the  collision.  In 
addition,  let  m  and  -e  be  the  mass  and  charge  of  the 
molecular  electron. 

During  the  collision,  the  center  of  mass  of  the  two- 
particle  system  moves  with  a  uniform  velocity  given  by 

vc ■  -  m  Vo  ( A—  1  ) 

M  +  m 

In  the  center  of  mass  (c.m.)  system  the  particles  move  in 
hyperbolic  orbits  (see  Fig.  A- 1 ) .  Let  vi  and  bi  be  the 
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initial  velocity  and  impact  parameter  of  M  relative  to  the 
c.m.  We  then  have 

vi  =  m  v0  (A-2) 

M  +  m 

and 

bi  =  m  b  ( A-3 ) 

M  +  m 

Let  ri  be  the  distance  of  M  from  the  c.m.  at  any  time.  Then 
the  force  acting  on  M  is  given  by 

F  =  e  E  m2  (A-4) 

(M  +  m)*  ri! 

and  is  directed  towards  the  c.m. 

If  9  is  the  angle  between  the  asymptotes  of  the 
hyperbola  then 

tan  (9/2)  =  c/a  (A-5) 

where  a  and  c  are  the  axes  of  the  hyperbola.  It  is  a 

property  of  hyperbolae  that  c  =  bi  .  In  addition,  if  J4  is 
the  acceleration  of  M  at  a  unit  distance  from  the  c.m., 

X  =  e  E  m2  ( A  —  6 ) 

M  (M  +  m ) * 


then  vi  1  =  ft/a ,  or 


a  = 


e  E  m2 


vi  2  M  (  m  +  M ) : 


(  A  —  7  ) 


Using  Eqns.  A-5  and  A-7  we  then  have 
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( A— 8 ) 


tan  (8/2)  =  bi  vi  »  M  (M  +  m )  * 

e  E  m* 


or,  with  Eqns .  A-2  and  A-3, 

tan  (8/2)  =  b  va »  m  M  (A-9) 

e  E  (M  +  m) 

After  the  collision,  the  component  of  the  velocity  of  M 
parallel  to  v0  is 


v»  =  M  va  -  m  Vo  cos  8  ( A—  1 0 ) 

M  +  m 


and  the  component  perpendicular  to  v0  is 

v„  =  m  va  sin  8  ( A— 11) 

M  +  m 

Therefore,  the  total  kinetic  energy  of  M  after  the  collision 
is  given  by 


Ek  =  1/2  M  vt  *  +  1/2  M  vB* 


or 


Ek  =  i  M  va*  (  M*  +  ra«  -  2  m  M  cos  8  ]  ( A- 1 2 ) 

2  (M  +  m)* 

The  energy  transfered  to  the  molecular  electron  is  then 
the  change  in  kinetic  energy  of  the  incident  particle: 


Q  =  1/2  M  v0J  -  Ek 
or 

Q  =  (  1  +  cos  8)  n  M1  v0  2  ( A- 1 3  > 

(M  +  m)» 
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Now 


1  +  cos  9=2  cos1 (  0/2 )  =  2  [1  +  tan2 { 0/2 ) ] * 1 

so  we  get,  using  Eqn.  A-9, 

Q  =  2  e2  E2  1  ( A- 14 ) 

m  Vo 2  b2  +  d2 

where 


d  =  e  E  ( M  +  m  1 
m  M  Vo 2 


Eqn.  A-14  thus  represents  the  transfer  of  energy  from  the 
incident  particle  to  the  molecular  electron. 

It  should  be  noted  at  this  point  that  the  above 
calculations  were  made  with  the  assumption  that  the  electrons 
in  the  molecule  are  free  to  move.  If  the  forces  holding  the 
molecular  electron  are  small  enough  that  the  period  for  a 
small  oscillation  of  the  electron  about  an  equilibrium  point 
is  large  compared  to  the  duration  of  the  collision,  then  the 
electron  can  be  considered  free.  Therefore,  to  be  able  to 
use  the  above  energy  transfer  calculation,  the  incident 
particle  must  be  sufficiently  fast.  A  90  V  electron  is 
moving  at  about  6  x  10®  cm/sec.  If  the  distance  over  which  a 
collision  occurs  is  on  the  order  of  one  Bohr  radius,  then  the 
collision  time  is  approximately  10*  17  sec  -  short  enough  that 
this  classical  approach  should  be  valid  for  the  present  work. 

Now  that  the  energy  transfer  in  the  collision  is  known, 
it  is  necessary  to  find  the  energy  distribution  of  the 
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molecular  electrons.  The  chance  of  a  transfer  of  energy 
between  Q  and  Q  +  dQ,  f(Q)dQ,  is  equal  to  the  chance  2Inb  db 
that  b  is  between  b  and  b  ♦  db,  n  being  the  number  density  of 
electrons  in  the  molecule.  Therefore,  we  have 

f(Q)  dQ  =  2  IT  n  b  db  (A-15) 

now  2  b  db  =  db1 ,  so 

f ( Q )  dQ  =  1  n  dbf.  dQ  (A- 16) 

dQ 


From  Eqn.  A-14  we  get 


db*  =  -  2  e2  E2  1  (A-17) 

dQ  a  Vo 2  Q2 

Putting  this  in  Eqn.  A-16  and  ignoring  the  minus  sign  we  have 

f ( Q )  dQ  =  2  1  n  e2  E2  1  dQ  (A- 18) 

m  v0 2  Q2 


Therefore,  for  an  incident  electron,  the  distribution 
for  the  molecular  electron  is 


f(Q)  =  2  1  n  e4 

m  Vo  2  Q2 


(  A-  1 9  ) 


If  we  regard  the  production  of  secondary  electrons  in  a 
MMIS  as  a  case  of  ionization,  we  are  then  saying  that  an 
electron  emitted  from  the  cathode  will  collide  with  an  Hj 
molecule  causing  ionization.  The  molecular  electrons  will 
then  have  a  distribution  given  by  Eqn.  A-19.  Since  the 
secondary  electrons  have  to  do  work  equal  to  the  ionization 
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energy,  Ei ,  to  escape  the  molecule,  the  distribution  of 
energy  for  the  secondary  electrons  will  be  given  by 


ft  boa(E) 


(E  ♦  Ei  )* 


( A-20 ) 


The  constant  A  can  be  found  by  requiring  normalization 
of  the  above  distribution: 


/  A  (E  +  Ei  )- 1  dE  =  1 
o 


( A-21 ) 


This  gives,  with  straight-forward  integration, 


A  =  Ei  (  1  +  Ei  /  V) 


( A-22 ) 


Thus,  the  Thomson  distribution  of  secondary  electrons  is 


ft  h o  ■  (  E )  =  Ei  M  t  Ei  /V) 
(E  +  Ei  ) 2 


(30) 
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APPENDIX  B  -  MOMENT  BQUATION  CODE 


In  this  appendix,  the  numerical  solution  of  the  moment 
equations  is  discussed.  After  this  a  listing  of  the  computer 
code  is  given.  A  brief  schematic  of  the  numerical  process  is 
given  in  Fig.  B-l. 

Calculation  of  Reaction  Rates 

Maxwell ian-averaged  reaction  rates  as  a  function  of 
average  electron  energy,  defined  by  Eqn.  4,  were  calculated 
once  at  the  beginning  of  the  study  and  stored  in  ASCII  files. 
These  files  are  then  referenced  each  time  the  main  computer 
code  is  run. 

In  order  to  calculate  the  reaction  rates,  the  (Ji  (E)v(E) 
data  was  fit  using  natural  cubic  splines  [30:122].  The 
integration  in  Eqn.  4  was  then  carried  out  over  the  range  of 
the  data,  for  average  electron  energies  ranging  between  .01 
eV  and  1000  eV .  The  integration  was  performed  using 
Simpson’s  adaptive  quadrature  [30:174-175]. 

The  Numerical  Solution  of  the  Moment  Equations 

The  main  computer  code,  MOMEQN ,  uses  the  above 
calculated  Maxwellian  reaction  rates  in  order  to  calculate 
the  time  evolution  of  the  number  densities  and  electron 
temperature  by  simultaneously  solving  the  moment  equations, 
represented  by  Eqns .  43-47  and  50. 

As  further  input  to  the  program,  MOMEQN  reads  two  ASCII 
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files  -  'process’,  and  ’ init*.  The  file  ’process’  contains  a 
line  of  data  for  each  plasma  reaction  that  waB  considered  n 
this  model.  The  first  element  of  each  line  is  a  character 
string  mnemonic  refering  to  the  reaction.  The  second  element 
is  either  a  1  or  a  0,  indicating  whether  (1)  or  not  (0)  that 
particular  reaction  is  active  during  the  current  run.  The 
third  element  is  a  real  number  equal  to  the  threshold  energy 
for  the  reaction.  The  file  ’init’  contains  all  of  the  other 
operating  data,  such  as  the  discharge  current,  voltage, 
pressure,  plasma  volume,  initial  values  of  the  densities, 
etc . 

After  all  of  the  above  data  is  entered  into  the  program, 
MOMEQN  calls  the  subroutine  DTAIL  to  calculate  the  high 
energy  tail  of  the  EEDF.  This  is  accomplished  by  first 
fitting  the  Oi  (E)v(E)  data  with  natural  cubic  splines.  Then 
ft  ii  i  (E)  is  calculated  using  Eqn .  32,  starting  at  E  =  V  -  Ei  . 
The  rest  of  the  distribution  is  calculated  for  successively 
lower  energies,  using  the  knowledge  of  what  the  distribution 
is  for  higher  energies. 

Once  this  distribution  is  calculated,  MOMEQN  calls  the 
subroutine  ELECTRONIC  in  order  to  calculate  the  reaction 
rates  averaged  over  the  tail  of  the  distribution,  which  will 
be  used  in  the  vibrational  density  equation.  These  reaction 
rates  are  determined  by  using  Simpson’s  adaptive  quadrature 
to  perform  the  integration  in  Eqn.  39.  The  particular 
reactions  for  which  rates  are  calculated  at  this  point  are 
the  electronic  excitations  and  H2  ionization. 
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Next,  MOMEQN  calls  the  subroutine  THOMSON,  which 
calculates  the  average  energy  of  the  source  (Eqn.  53)  and  the 
average  probability  of  no  wall  collisions  (Eqn.  59).  These 
integrals  are  again  calculated  using  Simpson’s  adaptive 
quadrature . 

At  this  point,  MOMEQN  is  ready  to  start  the  integration 
of  the  moment  equations.  This  is  accomplished  by  first 
setting  the  time  variable  to  0,  and  then  calling  the 
subroutine  RKF45.  This  subroutine  integrates  the  moment 
equations  from  the  previous  time  to  the  present  time  and 
outputs  the  answer.  The  time  variable  is  then  incremented 
and  RKF45  is  called  again.  This  process  continues  until  one 
of  two  conditions  is  met:  1)  the  number  densities  and 
electron  temperature  all  have  reached  a  steady  state,  or,  2) 
the  time  variable  has  reached  a  maximum,  user-specified, 
value.  Condition  1  is  numerically  determined  to  be  achieved 
if  the  values  of  all  the  number  densities  and  the  electron 
temperature  have  not  changed  by  more  than  a  set  factor  (10' 4  ) 
for  four  consecutive  time  steps.  If  this  condition  is  met, 
the  values  have  reached  a  steady  state,  as  far  as  the 
computer  code  is  concerned. 

The  subroutine  RKF45 

The  subroutine  RKF45  integrates  a  system  of  coupled 
first  order  differential  equations  by  using  two  Runge-Kutta- 
Fehlberg  integration  schemes.  One  is  of  order  4  and  the 
other  is  order  5.  The  two  results  obtained  are  compared  to 


get  an  estiiate  of  the  error  in  the  answer.  This  estimate  is 
then  used  for  step  size  control. 

The  Runge-Kutta-Fehlberg  algorithm  calculates  the  value 
of  the  various  plasma  parameters  by  using  the  following 
formula : 


6 

yj .  in  =  yj.«  +  sJfi  kj.i 

i  =  i 


( B- 1  ) 


where 

i  -  1 

kj  I  1  =  h  a  y’j(yjiB  +  2  Bin  kj  ,  ■  ,  tB  hn)  (B-2) 

■  =  i 

and  where  yj  ,  B  is  the  estimate  for  the  j-th  plasma  parameter 
(electron  density,  ion  densities,  or  electron  temperature)  at 
the  n-th  time  step,  y’  is  the  corresponding  derivative,  hB  is 
the  step-size  for  the  n-th  time  step,  tB  is  the  time,  and  cX  , 
Q,  and  %  are  coefficients  that  were  found  by  Fehlberg.  He 
found  one  set  of  s  that  makes  the  estimate  given  by  Eqn .  B- 
1  order  4  and  another  set  of  that  makes  it  order  5. 

The  majority  of  the  code  in  RKF45  is  used  to  determine 
the  optimum  step-size,  h,  for  the  current  time. 

Implementation 

The  program  MOMEQN  was  written  in  FORTRAN  77  and  ran  on 
a  VAX  11/785  under  Berkeley  4.3  BSD  UNIX.  Typical  run-times 
varied  between  10  minutes  and  1  hour,  depending  upon  the 
system  load.  It  is  quite  possible  that  many  improvements 
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could  be  made  on  this  code  to  make  it  run  more  efficiently. 
The  present  form  of  the  code  was  written  for  readability. 


( 

»* 


A  SELF-CONSISTENT  MOMENT  EQUATION  SOLVER  FOR  A 
MAGNETIC  MULTICUSP  HYDROGEN  DISCHARGE 


Written  by  David  E.  Bell,  2Lt,  USAF 


SYNOPSIS:  This  program  solves  the  time-dependent 

moment  equations  developed  in  Chapter  II  for  a  given 
set  of  initial  conditions.  It  outputs  the  time 
evolution  of  the  ionic  species  and  the  electron 
temperature . 


VARIABLE  LIST: 


c  on  -  array  indicating  whether  or  not  (1/0)  a  process  is  on 
c  th  -  array  containing  the  thresholds  for  the  processes  (eV) 
c  press  -  pressure  of  the  discharge  (Torr) 
c  volt  -  discharge  voltage  (V) 
c  curr  -  discharge  current  (A) 
c  dt  -  time  increment  for  output  (sec) 
c  y  -  array  whose  elements  represent: 
c  1  -  number  density  of  electrons  (cm~-3) 

c  2  -  number  density  of  negative  hydrogen  ions  (cm'-3) 

c  3  -  number  density  of  positive  H2  ions  (cm~-3) 

c  4  -  number  density  of  positive  H  ions  (cmA-3) 

c  5  -  average  electron  temperature  ( eV ) 

c  6  -  number  density  of  vibrationally  excited  H2  (cm'-3) 

c  nho  -  ratio  of  H  to  H2  in  plasma 
c  no  -  density  of  H2  (cm'-S) 
c  nh  -  density  of  H  (cra~-3) 
c  so  -  source  strength 
c  ens  -  average  energy  of  the  source 
c  ml  -  mass  of  an  electron  (kg) 
c  m2  -  mass  of  a  hydrogen  atom  (kg) 
c  t  -  time  variable  (sec) 
c  tout  -  next  time  of  output  (sec) 
c  tfinal  -  time  to  quit  program  (sec) 
c  temp  -  H2  temperature  (eV) 

c  relerr  -  relative  error  allowed  in  output 
c  abserr  -  absolute  error  allowed  in  output 
c  work  -  array  used  by  rkf45  for  various  things 
c  vol  -  plasma  volume  of  source  (cm~3) 

c  area  -  effective  wall  loss  area  of  plasma  source  (cm”2) 
c  vp  -  plasma  potential  (V) 

c  iwork  -  array  used  by  rkf45  for  various  things 
c  iflag  -  flag  used  by  rkf45  to  indicate  progress.  Values: 
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c  -2  -  single  successful  step  in  direction  of  tout 

c  2  -  integration  reached  tout.  Successful  return, 

c  3  -  relative  error  too  small,  has  been  increased 

c  4  -  more  than  3000  derivative  evaluations  needed 

c  5  -  solution  vanished  -  cannot  test  relative  error, 

c  absolute  error  must  be  nonzero  to  continue 

c  6  -  requested  accuracy  unobtainable,  error 

c  tolerance  increased 

c  7  -  RKF45  is  inefficient  for  this  problem 

c  8  -  invalid  input  parameters  : 

c  neqn.lt.O 

c  t=tout  and  iflag.ne.+l  or  -1 

c  relerr  or  abserr  .It.  0 

c  iflag.eq.O  or  .It. -2  or  .gt.8 

c  neqn  -  number  of  differential  equations  in  system  (6) 
c  file  -  character  array  with  file  names  for  the  processes 
c  sav  -  array  used  to  save  y( )  data  for  checking  convergence 
c  flag  -  logical  variable  used  in  testing  convergence 
c  ss  -  counter  used  in  testing  convergence 
c  length  -  average  distance  to  wall  (cm) 
c  xsect  -  sigma*vel  data  for  electronic  excitations 
c  nxsect  -  number  of  data  points  for  each  excitation  in  xsect 
c  a,b,c,d  -  cubic  spline  data  for  xsect 

c  tail  -  reaction  rate  coefficient  for  electronic  states 
c  averaged  over  tail  of  distribution 

c  rates  -  array  holding  Maxwellian  reaction  rate  data 
c  ftail  -  array  representing  the  tail  of  the  distribution 
c 

c ttt***t**t*ttttt*%tt%%t%tt**%t***t**tt*ttttttt*t**ttt*****tt 

external  diffeq 

real  on( 27 ) , th ( 27 ) , press , volt , curr ,dt , y ( 6 ) , nho , no ,nh , 
&  so, ens ,ml ,m2 , t , tout , tf inal , temp , relerr , abserr , 

&  work ( 39 ) , vol , rates (27,561,2), area , vp , sav ( 6 ) , 

&.  length  ,xsect(9,50,2)  ,  a  (  9 , 50  )  ,  b  ( 9 , 50  )  ,  c  (  9 , 50 )  , 

&  d ( 9 , 50 ) , tail ( 8 ) ,ftail(0: 100) 

integer  i work (5) ,iflag,neqn,ss,nxsect(9) 
logical  flag 
character*7  file(27) 

common  so,no,nh 
common /ener/ens , m 1 , m2 
comraon/onoff /on , f i le 
common/ freq/th ,rates,tail 
common /geom/ area , vol , vp 
common/ recom/y , temp , length 
common/ power /volt , curr 
coramon/s igv/xsect , nxsect 
common/ spl in/ a , b , c , d 
common/dist/f tai 1 

c 

c  Input  rate  coefficient  data  from  files 


call  process ( on , th , f i le , rates , xsect , nxsect ) 

Input  initial  operating  conditions 

call  ini t( press , volt , cur r , tf inal ,dt #y(l),y(2),y(3), 
&  y ( 4 ) , nho , temp , vol , area , length , vp , re 1 err , 

&  abserr ) 

no =8 . 3216el4*press/temp/( 1 . fnho) 
nh=nho*no 
sav( 1 ) =y( 1 ) 
y( 1 )  =  1 .  ell 


Calculate  tail  of  EEDF 
call  dtail(ftail) 

Calculate  rate  coefficients  for  electronic  excitations 
call  electronic 

Calculate  average  energy  of  source  and  source  strength 

call  thomson ( ens , volt , so ) 
y (  1 )  =  sav( 1 ) 

so=so*curr/l . 602 2e- 19/ vol* volt/ ( 2 .*15.427) 

y ( 5 ) =ens 

ml  =  9 . 1095e-3 1 

m2=2.*1.6726e-27 

Set-up  output  files 

call  files 


c  Start  integration  of  equations 
c 

t  =  0 . 
i  f  lag=  1 
tout  =  t 
neqn  =  6 

10  call  rkf 45 ( diffeq , neqn , y, t , tout , relerr , abserr , i flag 

&  work , iwork ) 

call  output(t,y) 

goto  (80,20,30,40,50,60,70,80) , if lag 
c 

c  Successful  step  in  integration 
c 

20  tout=t+dt 
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c  Check  for  convergence  of  answer 
c 

f lag= . true . 
do  21  i= 1 , 6 
if (y(i) .eq.O. )then 

f lag= ( abs (sav(i)-y{i)).lt. abserr/ 10.). and .flag 
else 

flag=(abs(sav(i)-y(i) )/y(i) .lt.relerr/10. ) . and .flag 
endif 

21  continue 

i f ( f lag ) then 
ss  =  ss+  1 
else 
ss  =  0 
endif 

do  2  3  i  =  1 , 5 
23  sav(i)=y<i) 


c  Four  consecutive  times  of  little  change  in  output 
c  constitutes  convergence 


if ( ss . eq . 4 ) goto  22 
if ( t . It . tf inal )goto  10 


c  Completed  successful  run  of  program  -  output  summary 
c 

22  call  endout(y) 

stop 


c  Errors  in  RKF  routine  handled  here 


write{6,31)t,relerr, abserr 

format! 1  Tolerances  reset  at  time  ’ ,lpe9.3/ 

’  Relative  error  =  ’,lpe8.2/ 

’  Abserr  =  ’ , lpe8 . 2 ) 

goto  10 
wri te ( 6 , 4  1  )  t 

format ( lx , ’Over  3000  function  evaluations  at  time 
lpe9 .3) 

goto  10 
abserr: 1 . e-9 

write ( 6 , 31 ) t , relerr , abserr 
goto  10 

relerr=10. *relerr 
write(6,31 )t,relerr,abserr 
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if lag  =  2 
goto  10 


70 

71 


& 


write ( 6 , 7 1 ) t 

f ormat ( lx , ’ RKF4 5  is  inefficient  for  this  problem.’/ 
’Time  =  ’ , lpe9 . 3 ) 

if lag=2 
goto  10 


80 

81 


write(6,81 )t 

format ( lx ,  ' Improper  call  at  time  ’,lpe9.3) 
stop 


end 


******** 

c 

This  subroutine  opens  the  file  ’process’  and  reads 
it  to  determine  what  processes  are  active  for  the 
current  run  (on),  the  threshold  energy  for  the 
processes  (th),  and  the  names  of  the  files  that 
contain  the  rate  look-up  data  (file).  The  rate  data 
is  then  entered  into  the  array  rates()  which  is  used 
by  the  function  rate()  to  determine  the  rate 
coefficient  for  a  given  energy.  This  subroutine  also 
enters  the  cross-sect ion*velocity  data  as  a  function 
of  energy  for  the  electronic  excitations. 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


The  file  ’process’  is  made  up  of  27  records  of  the 
form : 


file 


1  .  ( on/of  f ) 


12.0  (threshold) 


c************************************************************ 


subroutine  process (on, th, file, rates, xsect.nx sect) 
real  on(27) , th ( 27 ) ,rates(27,561,2) ,d2,d3,xsect(9,50,2) 
integer  nxsect(9) 
character*?  file(27) 
character*40  dl 


open(unit=15,file='process' ,status=’old’ ) 


do  1  i  =  1  , 27 

read( 1 5 , * ) f i le ( i ) ,on( i ) ,th(i) 


c  lose ( un i t= 1 5 ) 
do  5  i  =  1  ,27 


c 

c 

c 


the  extension  .rat  is  used  to  indicate  average  rate  data 
open(unit=15,file  =  file(i)//’  .rat’  ,status=’old’  ) 
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fi 


a 


01 


a 


a 


V! 


4 


4 


'A 


first  five  lines  of  files  are  descriptive  labels  used  for 
graphing  purposes  only 
sixth  line  contains  the  number  of  data  pairs 
seventh  line  contains  xmin,xmax 
eigth  line  contains  ymin,ymax 

do  2  j=l , 5 
read!  15 , * )dl 
read( 15 , 3 )n 
format (id) 
read( 15,*)d2,d3 
read{ 1 5 , * )d2 , d3 

do  4  j  =  l,n 
read( 15 , * )d2 ,d3 
rates ( i , j , 1 ) =d2 
rates ( i , j , 2 ) =d3 
continue 

close ( unit= 1 5 ) 

continue 

do  6  i=l ,9 


the  extension  .vel  is  used  to  denote  cross¬ 
sect  ion*velocity  data 

open ( unit=15,file=file(i)//’ .vel’ ,status=’old’ ) 

do  7  j  =  1 , 5 

read( 1 5 , * )dl 

read ( 1 5 , 3 ) nxsect ( i  ) 

read) 15 , * )d2 , d3 

read( 15 , * )d2 ,d3 

do  8  j= 1 , nxsect ( i  ) 

read( 15,*)xsect(i,j,l),xsect(i,j,2) 

close! uni  t  = 1 5  ) 

conti nue 

return 

end 


This  subroutine  opens  the  file  ’init’  and  reads  in 
the  initial  operating  conditions  of  the  discharge. 


»•  t*> 


1 
» 
t 

t 

subroutine  ini t ( press ,  volt , curr , tf inal ,dt , ne , nm , np,  j 

nhp , nho , temp , vol , area , length , vp ,  j 

relerr , abserr )  * 

real  press , volt ,curr , tf inal , dt , ne , nm , np , nhp, nho , temp ,  j 

&  vol , relerr , abserr , area , vp , length  j 

open ( unit=15 , f ile= * init ’ , status= ’ old’ )  \ 

I 

read( 15 , * (press  ! 

read (15,*) volt 
read( 15 , * (curr 

read( 15 , * )tf inal  ; 

read( 15 , * )dt  ! 

read ( 1 5 , * ) ne  ! 

read( 15 , * )nm  J 

read ( 1 5 , * ) np  ] 

read (15,* (nhp 
read( 15 , * (nho 

read( 15 , * ) temp  ; 

read( 15 , * (vol 

read( 15 , * (area  I 

read( 15 ,*) length 

read( 15 , * ) vp 

read( 15 , * (relerr 

read( 15 ,* (abserr 

close ( unit=15 ( 

return 

end 

c ************************************************************ 

c 

c  This  subroutine  sets  up  the  output  files  for  the 

c  number  densities  of  tne  ionic  species  and  the 

c  electrons  and  the  electron  temperature.  They  are 

c  set  up  in  the  correct  format  to  be  used  by  'graph', 

c  The  program  'graph'  takes  a  data  file  and  produces  a 

c  file  that  :1s  in  the  correct  format  to  be  plotted 

c  either  on  a  Benson  printer  or  on  a  Laser  printer, 

c  The  first  5  lines  of  the  data  file  must  be  blank  or 

c  contain  labels  for  the  plot.  The  next  line  must 

c  contain  the  number  of  data  points  in  i3  format.  The 

c  next  two  lines  are  the  minimum  and  maximum  values  for 

c  the  x  range  and  the  y  range.  The  rest  of  the  lines 

c  contain  tne  data  points. 

c  Another  file  ('all')  is  set  up  for  output  in  a  more 

c  readable  format, 

c 

c*tt*t**ttttt*tttttt*tttttt*tt****t***%**t*ttt*****ttt*tt*t%t 

subroutine  files 
character*!  line(80) 
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open (unit=10,file=’ne’ , status = ’ new’ ) 
open ( unit= 1 1 , f ile= ’ nra ’ , status= 'new* ) 
open( unit=12,file=’np’ , status = ’ new ’ ) 
open(unit=13 , f ile= ’nhp* , status= 'new* ) 
open ( uni t= 14, file= ’en’ ,status= ’new* ) 
open (unit=16,file='nv’ , statuB= ’ new’ ) 
open( unit=20 , f ile= ’ all ’ , status= ’new’ ) 

write ( 10 ,*)’ Electron  density’ 

write ( 1 1 ,*)’ H-  density’ 

write( 12, * ) ’H2+  density’ 

wri te( 13 , * ) ’ H+  density’ 

write ( 14 , * ) ’ Electron  temperature ’ 

write( 16 ,*) ’Vibrational  density* 

do  1  i  =  10 , 13 
write ( i , * ) ’ Time ’ 
write( i , * ) '  ( sec  )  ’ 
write ( i , * ) ’ Density ’ 
write(i,*)’  (cm~-3)’ 
continue 

write( 18 , * ) ’Time ’ 
write  <  16 , * )  ’  ( sec  )  ’ 

write ( 16 , * ) ’ Density’ 
write(16,*)’  (cm,'-3)’ 

write( 14 , * ) ’Time ' 
write ( 14 , * )  ’  <  sec )  ’ 
wri te ( 14 , * ) ’ Temperature ’ 
wri te ( 14 , *  )  ’  ( eV ) ’ 

do  3  i  =  10 , 14 
write ( i , 2 ) 999 
format< i3 ) 
write ( i , * ) 0 . , 0 . 
write ( i , * ) 0  .  , 0  . 
continue 

write ( 20 , 4 ) 

format ( ’TIME’ , 9x , ’ Ne ’ , 1  lx , ’ N- ’ , 1 lx , ’ N+ ’ , llx , ’Nh+’  , 
lOx , ’Te’ ) 
do  6  i= 1 , 80 
1 ine ( i ) = ’ = ’ 
continue 

write(20,5)(line(i) , i  =  1 , 8  0 ) 
format ( 80a ) 

return 

end 


***** 


This  subroutine  calculates  the  rate  coefficients  for 
the  electronic  excitations  based  on  the  modified 
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c 

c *** t t**t***t*t %%***** ttttttttttt *********** ttttttttt******** 

subroutine  electronic 

real  xsect(9,50,2) ,tol,a(9,50) , b( 9 , 50 ) , c( 9 , 50 ) , d( 9 , 50)  , 
k  rates (27, 56 1,2) , tl ( 200 ) , a2 ( 200 ) , h2 ( 200 ) ,fa(200) , 

k  fc(200) ,fb(200) , s ( 200 ) , k2 ( 200 ) ,fd,fe,sl,s2,v(8) , 

k  th ( 27 ) , ion, volt, curr , vol , area , vp , tail ( 8 ) 

integer  nxsect(9) 

common/ spl in/a , b , c , d 
common/ sigv/xsect , nxsect 
common/ f  req/th , rates , tai 1 
common/geom/area , vol , vp 
common/ power /volt , curr 

tol  =  l  .e-3 
ion=15.427 

do  1  k=  1 , 8 

tai 1 ( k ) =0 . 

i  =  1 

tl ( i ) = 10 . *tol 

a2 ( i ) =xsect (k , 1,1) 

h2 ( i ) = ( xsect ( 1 , nxsect (k) ,l)-xseot(k,l ,1) )/2. 
fa(i)=felect(k,xsect(k, 1 , 1) ) 
fc(i)=felect(k,xsect(k,l, l)+h2(i) ) 
fb( i)=felect(k,xsect(k,nxsect(k) ,  1)  ) 
s(i)=h2(i)*(fa(i)+4.*fc(i)+fb(i) )/3. 
k2 ( i )  =  1 . 

5  i f ( i . le . 0  )  goto  4 

if(tl(i).lt.i.e-6)tl(i)=l.e-6 
f d= f elect (k,a2( i )+h2( i)/2. ) 
fe=felect(k,a2(i)+3.*h2(i)/2. ) 
sl=h2( i ) * ( fa( i )+4 . *fd+fc( i ) )/6 . 
s2  =  h2 (  i  )  * ( fc ( i ) +4 . *fe  +  fb( i ) )/6  . 
v(l)=a2(i) 
v  (  2  )  =  f  a  (  i  ) 
v( 3 )  =  fc( i  ) 
v  (  4  )  =  f  b  (  i  ) 
v  (  5  )  =h2 ( i ) 
v  (  6  )  =  1 1  (  i  ) 

V ( 7 ) =8 ( i  ) 
v ( 8 ) =k2 ( i  ) 

i  =  i  -  1 

i f ( abs (sl+s2-v(7) ) . 1 t . v ( 6 ) )then 
tail(k)=tail(k)+sl+s2 
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else 

if ( v( 8 ) .  ge . 200 ) then 

write ( 6 ,*)’ Integration  failed  in  ELECTRONIC* 
stop 
else 
i  =  i  +  l 

a2 ( i ) =v( 1 ) +v( 5 ) 
fa( i ) =v( 3  ) 
f c ( i ) =f e 
fb( i ) =v( 4  ) 
h2(i)=v(5)/2. 
tl(i)=v(6)/2. 
s( i ) =s2 
k2 ( i ) =v( 8 ) +1 . 

i  =  i  +  l 

a2 ( i ) =v{ 1 ) 
fa( i ) =v( 2 ) 
fc( i ) =fd 
fb( i ) =v( 3  ) 
h2 ( i ) =h2 ( i  —  1 ) 
tl(i)=tl(i-l) 
s ( i ) =sl 
k2 ( i ) =k2 ( i- 1 ) 
endif 
endif 

goto  5 

4  continue 

3  continue 

1  continue 


return 

end 


c************************************************************ 


c 

c 

c 

c 


This  subroutine  performs  natural  cubic  splines  on  the 
cross-section  data  for  the  electronic  excitations. 


20 

10 


subroutine  spline (data, nun ,a,b,c,d) 
real  data(9,50,2) , a ( 9 , 50 ) , b ( 9 , 50 ) , c ( 9 , 50 ) , d ( 9 , 50 ) , 
&  h ( 50 ) .alpha (50) , 1 ( 50 ) , mu ( 50 ) , z ( 50 ) 

integer  num(9) 

do  10  k=l ,9 
do  20  j= 1 , num( k ) 
a( k  ,  j ) =data ( k , j , 2 ) 
continue 
continue 
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do  1  k= 1 , 9 


do  2  i=l, num ( k ) - 1 

2  h( i ) =data(k , i  + 1 , 1 ) -data (k ,  i  ,  1 ) 
do  3  i=2,num(k)-l 

alpha { i)=3.*(a(k,i+l)*h(i-l)-a(k,i)*(data(k>i+l,l) 

&  -data(k,i-l, 1) )+a(k,i-l)*h(i) )/(h(i-l)*h(i) ) 

3  continue 

1(0)=1. 
mu ( 0 ) =0 . 
z ( 0 ) =0 . 

do  4  i=2,num(k)-l 

l(i)=2.*(data(k,i+l,l)-data(k,i-l,l))-h(i-l)*mu(i-l) 

mu(i)=h(i)/l(i) 

z  (  i  )  = ( alpha ( i)-h(i-l)*z(i-l) )/l(i) 

4  continue 

1 ( num( k ) ) = 1 . 
z(num(k) ) =0 . 
c  (  k , num ( k ) ) =0 . 

do  5  j=num(k ) -1 , 1 , -1 

c(k , j )=z( j ) -mu ( j ) *c(k , j+1 ) 

b(k,j)=(a(k,j+l)-a(k,j))/h(j)-h(j)*(c(k,j+l) 

4  +2.*c(k, j) )/3. 

d(k,j)=(c(k,j+l)-c(k,j) )/(3.*h(j) ) 

5  continue 

1  continue 


return 

end 


c 

c  This  function  represents  the  distribution  function 

c  for  the  tail  multiplied  by  sigraa*velocity . 

c 


real  function  felect(i,e) 
real  e , i on , f ta i 1 ( 0 : 1 00 } , di s t r ib 

common/dist/f tai 1 

ion:  15 . 427 

distrib=(e-float( int(e) ) )  * 

4.  (ftail(int(e)  +  l)-ftail(int(e)))  +  ftail(int(e)> 
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felect=sv( i,e)*distrib 

return 

end 

c t*tttttttttt*tttt*tttztttttttttttttttt*tttttt*****tttttttttt 

c 

c  This  function  uses  the  spline  data  calculated  by 

c  spline  to  calculate  sigma*vel  as  a  function  of  energy, 

c  This  routine  is  used  by  fcoli  and  felect. 


real  function  sv(i,en) 

real  en,xsect(9,50,2) ,den,a(9,50) , b( 9 , 50 ) , c ( 9 , 50 ) , 

&  d{ 9 , 50 ) 

integer  nxsect(9) 

common /s igv/xsect,nxsect 
common/ s pi in/a , b , c , d 

if((en.lt.xsect(i,l,l)) 

&.  . or . ( en . gt . xsect ( i , nxsect ( i ) , 1 ) ) ) then 

sv  =  0  . 
return 
endi  f 

do  1  j  =  1  ,  nxsect ( i ) -  1 

if(en.ge.xsect(i,j,l) . and . en . It . xsect ( i , j+ 1 , 1 ) ) goto  2 

1  continue 

j  =  nxsect ( i ) 

2  den=en-xsect ( i , j , 1 ) 

sv  =  a  (  i  ,  j  )  +b  (  i  ,  j  )  *den  +  c  (  i  ,  j  )  *den  *  *  2  +  d  (  i  ,  j  )  *  d«'  r  *• 

return 

end 

. . .  <  >  •  i  •  <  *  • 

c 

c  This  function  calculates  *  h»-  ’  '  •* 

c  rate  as  a  funet  ion  of  .* n •  i  *\ 

c 

etmtitutiitutmt  mum  m  ii  t  i  m  ■  i  >  i  i  o  .  > 

real  function  f  ■  .  .  ■  ■ 

real  en.no.s  ■  . nh 

i ■  i immi  in  s«  • .  n-  .  ?.► 
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fcoll=fcoll+sv(i,en) 

continue 


fcoll=no*fcoll 

return 

end 

c************************************************************ 

c 

c  This  function  calculates  the  rate  of  loss  of  primary 

c  electrons  to  the  walls  of  source.  Formula  is  from 

c  Bretagne,  et  al  [10:813]. 

c 

real  function  fwall(en) 
real  en , area , vol , no , vp, so , nh 

common  so,no,nh 
common/geom/area , vol , vp 

i f ( en . It . vp ) then 
fwall=0 . 
else 

fwall=6 . e7*area/vol*sqrt (en)»(l.-vp/en)/4. 
endif 

return 

end 

c************************************************************ 

c 

c  This  subroutine  calculates  the  tail  of  the 

c  distribution, 

c 

subroutine  dtail(ftail) 

real  ftail(0:100) ,volt,curr,th(27) , rates(27,561,2), 

&  tail(8),area,vol,vp,ion,a(9,50),b(9,50),c(9,50), 

Sl  d(9,50>,xsect(9,50,2),x,y,xp,yp,no,nh,so 

integer  nxsect ( 9 ) , f lag 

common  so,no,nh 
common /power /vol t , curr 
common/ freq/th, rates, tail 
common/geom/area , vo 1 , vp 
common/ s pi in/a ,b,c,d 
common/ s igv/xsect , nxsect 

call  spli ne(xsect, nxsect, a, b,c,d) 
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ion=15 . 427 
f lag=0 


do  1  i=100,0,-l 


en=f lost ( i ) 


if ( en. gt . volt-ion ) then 
f tail ( i ) =0 . 
goto  1 
endif 


sum=0 . 

do  2  j  =  1 , 9 

k=i+int ( th{ j ) + . 5 ) 

if (k.gt. 100)k=100 

sum=sum+f tai 1 (k)*sv(j,en+th(j) ) 

continue 

sum=no*sum 


x=fcoll(en) 

y=fwall(en) 

if ( ( f lag . eq. 0 ) . and. ( x+y .eq . 0 . ) ) then 
xp=fcoll ( en+ 1 . ) 
yp=fwall ( en+1 . ) 
f lag=  1 
endif 

if ( flag . eq. 1 ) then 
x=xp 

y=yp 

endif 

ftail(i)=(ion*l.+ion/volt)*curr/1.6022e-19/vol*volt 
/(2.*ion)/(en+ion) *  *2  +  sum ) / ( x  +  y ) 


continue 


return 

end 


This  subroutine  outputs  the  current  data  to  the 
output  files. 


(;******«**************«*********** *************************** 


subroutine  output(t,y) 
real  t , y ( 6 ) , telec 


write ( 10 , * ) t , y ( 1 ) 
write! 1 1 , * ) t , y ( 2 ) 
write! 12,*)t,y(3) 
write! 13 , * ) t , y( 4 ) 
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telec=2 .  *y ( 5 ) /3 . 
write( 14 , * ) t , telec 

write (20, l)t,y(l),y(2),y(3),y(4), telec 
1  format(6( lpe9.3,4x)  ) 

return 

end 

c 

c  This  subroutine  represents  the  system  of  differential 

c  equations  to  be  solved, 

c 

subroutine  dif f eq( t , y , yp ) 

real  t,y(6) ,yp(6) ,so,no,nh,ens,ml,m2,r(27) , th ( 27 ) 

real  rates ( 27 , 56 1 , 2 ) , tail ( 8 ) 

common  so,no,nh 

common/ener/ens , ml , m2 

common/ f  req/ th , rates , tai 1 

do  1  i =  1 , 2  7 

1  r( i ) =rate( i , y( 5 ) ) 

if ( y ( 5 ) .ge . 1000 . ) write( 6 , * ) ' Energy  >  1000  --  *,y(5) 
c 

c  electron  continuity 
c 

yp(l)  =  so  +  r(9)*y(l)*no  +  r (  19 ) *y ( 1 ) *no 
4  +  r( 13)*y( 1 )*nh  -  r< 18) *y( 1 ) *y(6 ) 

4  -  r ( 26 ) *y ( 1 ) * y ( 3 )  -  r( 15 ) *y( 1 ) *y( 3 ) 

4  -  r( 16 ) *y( 1 ) *y ( 4 )  +  r ( 22 ) *y ( 1 ) *y ( 2 ) 

4  +  r( 24 ) *y( 2 ) *nh  -  r ( 1 2 ) *y ( 1 ) *y ( 3 )  ♦  r ( 2 3 ) *nh*y ( 2 ) 

c 

c  negative  ion  continuity 
c 

yp(2)  =  r( 18) ty( 1 ) ty( 6)  ♦  r ( 25  )  *y ( 1 ) *no 
4  +  r(26)*y( 1 ) * y ( 3 )  -  r ( 2 1 ) *y ( 4 ) * y ( 2 ) 

4  -  r(22)<y( 1 ) *  y ( 2 )  -  r ( 24  )  *nh *y ( 2 ) 

4  -  r( 23 ) *nh*y( 2 ) 

c 

c  H2+  continuity 
c 

yp(3)  =  so  +  r(9)*y(l)*no  -  r (  1 2 ) *y ( 1 ) *y ( 3 ) 

4  -  r( 26 ) *y( 1 ) * y ( 3 )  -  r( 15 ) *y ( 1 ) *y( 3  ) 

c 

c  H+  continuity 
c 
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yp(4)  =  r{ 19 ) *y ( 1 ) *no  +  r ( 13 ) *y ( 1 ) *nh  +  r ( 25 ) *y ( 1 ) *no 
+  r ( 26 ) *y( 1 ) *y ( 3 )  -  r( 16 ) *y( 4 ) *y( 1 ) 

-  r(21)*y(4)*y(2) 


c  energy  equation 


yp ( 5 )  =  -  r( 1 ) *no*th( 1 )  -  r ( 2 ) *no*th ( 2 )  -  r ( 3 ) *no*th ( 3 ) 

-  r( 4 ) *no*th( 4 )  -  r ( 5  )  *no*th( 5 )  -  r ( 6  )  *no* th ( 6  ) 

-  r( 7 ) *no*th( 7 )  -  r ( 8 ) *no*th ( 8 ) 

-  r( 10 ) *no*th( 10 )  -  r ( 1 1 ) *no*th( 1 1 ) 

-  r ( 14 ) *no*th ( 14 )  -  r ( 9  )  *no* th ( 9 ) 

-  r< 13 ) *nh*th( 13 )  +  ens*so/y(l) 

-  r( 20 ) *2 . *ml/m2*y ( 5 ) *no  -  y( 5 ) /y ( 1 ) *yp( 1 ) 

yp( 5  )  =  yp( 5 )  -  r{12)*y(3)*y(5)  -  r { 15 ) *y ( 3 ) *y< 5 ) 

-  r(16)*y(4)*y(5)  -  r (  18 ) *y <6 ) *y < 5 ) 

-  r( 22 ) *y( 2 ) *y ( 5  )  -  r ( 23 ) *y ( 2 ) *nh/y ( 1 ) *y ( 5 ) 

-  r(24)*y(2)*nh/y( l)*y(5) 


vibrational  population  continuity 

yp(6)  =  (tail(l)  +  tail(2)  +  tail(3)  +  tail(4) 
k  +  tail(5)  +  tail(6)  +  tail(7)  +  tail(8))*no 

k  -  r( 18 ) *y( 1 ) *y ( 6 )  -  r ( 1 7 ) *y ( 6 ) / 10 . 

&  -  r ( 27 ) *y { 6 ) *nh 

return 

end 


c 

c  Fehlberg  Fourth-Fifth  Order  Runge-Kutta  Method 

c 

c  Written  by  H.A.  Watts  and  L.F.  Shampine 

c  Sandia  Laboratories 

c  Albuquerque,  NM 

c  As  published  in 

c  Forsythe,  G.E.,  Malcolm,  M.A.,  k  Moler,  C.B., 

c  COMPUTER  METHODS  FOR  MATHEMATICAL 

c  COMPUTATIONS,  Prentice-Hall,  NJ ,  1977, 

c  pp. 135-147. 

c 

c  This  subroutine,  in  conjunction  with  RKFS  and  FEHL, 

c  performs  a  fourth  order  and  a  fifth  order  Runge-Kutta 

c  routine  on  the  system  of  differential  equations.  It 

c  uses  the  output  from  these  two  calculations  to 

c  estimate  the  error  involved  in  the  answer  it  gives, 

c  If  the  error  is  too  great,  it  reduces  the  time  step 

c  and  tries  again.  The  process  continues  until  an 

c  answer  with  tolerable  error  is  returned  or  the  code 

c  decides  that  too  much  work  is  being  done.  For  a  more 

c  detailed  explanation  of  the  routine  see  the  above 


or 
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c  book. 

c  It  should  be  noted  that  the  routine  is  incorporated 

c  exactly  as  it  was  given  in  Forsythe, et  al .  As  such 

c  it  is  very  general.  It  is  possible  to  specialize  the 

c  routine  for  use  solely  by  this  code,  but  that  was 

c  deemed  unnecessary  and  undesirable  for  future 

c  expansion  of  the  model, 

c 

cttttttttttttttttttttttttttttt ******************************* 

subroutine  rkf45(f,neqn,y,t,tout,relerr,abserr,iflag, 

&  work,iwork) 

integer  neqn , i f lag , iwork ( 5 ) 

real  y ( neqn ) ,t,tout,relerr, abserr , work ( 3+6*neqn ) 
external  f 

integer  k 1 , k2 , k3 , k4 , k5 , k6 , k lm 

klm=neqn+ 1 

k 1 =k lm+ 1 

k2  =  k 1 tneqn 

k3=k2+neqn 

k4=k3+neqn 

k5=k4+neqn 

k6=k5+neqn 

call  rkf s  (  f  ,  neqn ,y,t,tout,relerr, abserr , i f lag , work (  1  )  , 

&  work(klm) ,work(kl ) ,work(k2) ,work(k3) ,work(k4 ) , 

&  work ( k5 ), work ( k6 ), work ( k6+ 1 ), iwork ( 1 ), iwork ( 2 ) , 

&  iwork ( 3 ) , iwork ( 4 ) , iwork ( 5 ) ) 

return 

end 

subroutine  rk fs(f, neqn, y,t, tout, re lerr, abserr , if lag , 

&  yp,h,fl,f2,f3,f4,f5,savre, savae , nf e , kop , 

&  int , j flag , kf lag ) 

logical  hfai Id , output 

integer  neqn, i f 1 ag , nf e , kop , int,jflag,kflag 

real  y(neqn) ,t,tout,relerr, abserr , h , yp( neqn ) , f 1 ( neqn ) , 

&  f 2 ( neqn )  ,  f  3  ( neqn ) , f  4 ( neqn ) , f 5 ( neqn ) , savre , savae 

external  f 

real  a,ae,dt,ee,eeoet,esttol,et,hmin,remin,rer,s, 

&■  scale, tol,toln,u26,epspl,eps,  ypk 

integer  k , maxnfe , raf lag 
real  amaxl,aminl 
data  remin/1 . e- 12/ 
data  maxnfe/ 3000/ 

i  f (neqn . It . 1 )goto  10 

if( ( re lerr. It. 0.  )  .or. (abserr. It. 0.  ) )goto  10 
mflag=iabs( if  lag) 

if( (mflag.eq.0) .or. (mflag.gt.8) )goto  10 
if(mflag.ne. 1 )goto  20 


Page  91 


•V 

.v 

.•» 


i 

I 

*  • 

£ 


VO 


■Si 


ssv 


5 


2? 


eps  =  1 . 
eps  =  epa/2 . 
epapl =eps+ 1 . 
i f ( epspl . gt . 1 . ) goto  5 
u26=26.*eps 
goto  50 

10  if lag=8 

return 

20  if ( ( t . eq. tout ) . and . ( kf lag . ne . 3 ) )goto  10 

if ( mf lag . ne . 2 ) goto  25 

if( (kflag.eq.3) .or. (int.eq.O) )goto  45 
if ( kf lag . eq . 4 ) goto  40 

if ( (kflag.eq.5) .and. ( abserr . eq . 0 . ) )goto  30 
if ( ( kf lag . eq. 6 ) . and . ( re 1 err . le . savre ) 

&  . and .( abserr . le . savae )) goto  30 

goto  50 

25  if ( if lag . eq. 3 )goto  45 

if ( if lag . eq. 4 ) goto  40 

i f ( ( i f lag . eq . 5 ) . and . ( abserr . gt . 0 . ) ) goto  45 

30  stop 

40  nfe=0 

if ( mf lag . eq . 2 )goto  50 

45  iflag=jflag 

i f ( kf lag . eq . 3 ) mf lag= iabs  (  i f lag ) 

50  jflag=iflag 

kf lag=0 

savre=relerr 

savae=abserr 

rer=2. *eps+remin 
i f ( relerr . ge . rer ) goto  55 

relerr=rer 
i f lag=3 
kf lag=3 
return 

55  dt=tout-t 

i f ( mf lag . eq . 1 ) goto  60 
i f { int . eq . 0 ) goto  65 
goto  80 

60  int=0 

kop  =  0 
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call  f(a,y,yp) 
nf  e=  1 

if { t . ne . tout ) goto  65 

if lag=2 

return 

int=  1 
h=abs (dt ) 
toln=0 . 

do  70  k=l,neqn 
tol=relerr*abs ( y ( k ) )+abserr 
if (tol.le.O. ) goto  70 
toln=tol 
ypk=abs ( yp( k ) ) 

if ( ypk*h**5 . gt . tol )h= ( tol/ypk ) ** .2 
continue 

if ( toln. le.O. )h=0. 

h=amaxl ( h , u26*amaxl ( abs ( t ) , abs ( dt ) ) ) 
jf lag=isign( 2 , if lag ) 

h=sign(h,dt ) 

if ( aba ( h ) . ge . 2 . *abs ( dt ) ) kop=kop+ 1 
if (kop.ne. 100)goto  85 

kop=0 
if lag=7 
return 

if ( abs ( dt ) . gt . u26*abs ( t )  )goto  95 

do  90  k= 1 , neqn 

y(k)=y(k)+dt*yp(k) 

a=tout 

call  f ( a , y , yp ) 
nfe=nfe+ 1 
goto  300 

output: . false . 

scale=2./relerr 

ae=scale*abserr 

hf ai ld= . false . 

hmin  =  u26*abs  (  t ) 

dt: tout- t 

i f ( abs (dt ) . ge . 2 . *abs ( h ) )goto  200 
i f ( abs ( dt ). gt . aba ( h )) goto  150 

output= . true . 
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h=dt 

goto  200 
150  h=.5*dt 

200  if { nf e . le . maxnfe ) goto  220 

if lag=4 
kf lag=4 
return 

220  call  fehl(f , neqn , y , t , h, yp , f 1 , f 2 , f 3 , f 4 , f 5 , f 1 ) 
nfe=nf e+5 

eeoet=0 . 
do  250  k= 1 , neqn 
et=abs(y(k) )+abs( f 1 (k) )+ae 
if ( et . gt . 0 . ) goto  240 

i f lag=5 
return 

240  ee=abs( ( -2090 . *yp{ k ) + ( 21970 . *f 3 ( k ) - 1 5048 . *f 4 ( k ) ) )+ 
&  (22528. *f2(k)-27360.*f5(k) ) ) 

250  eeoet=amaxl ( eeoet , ee/et ) 

esttol=abs(h ) *eeoet*scale/752400 . 

if(esttol.le.l.) goto  260 

hfai ld= . true . 
output= .false . 
s= .  1 

if ( esttol . It . 59049 . ) s= . 9/esttol** . 2 
h  =  s*h 

i f ( abs ( h ) . gt . hmin ) goto  200 

i f lag =6 
kf lag=6 
return 

260  t=t+h 

do  270  k= 1 , neqn 
270  y ( k ) =  f 1 ( k ) 

a  =  t 

call  f ( a , y , yp ) 
nfe=nfe+l 

s  =  5  . 

if(esttol.gt. 1. 88956 8e-4)s=. 9 /esttol**. 2 

if(hfaild)s  =  aminl ( s,  1  .  ) 

h=3 ign ( amax 1 ( s*aba ( h ) , hmin ) , h ) 

i  f  (  output ) goto  300 
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if ( if lag.gt.Oigoto  100 


if lag=-2 
return 

t=tout 
i f lag=2 
return 


subroutine  f ehl < f , neqn , y , t , h , yp , f 1 , f 2  ,  f 3  ,  f 4  ,  f 5  ,  s  ) 
integer  neqn 

real  y ( neqn  > , t , h , yp{ neqn ) , f 1 ( neqn ) , f 2 ( neqn ) , f 3 ( neqn ) , 
f 4 ( neqn ) , f 5 ( neqn ) , s ( neqn ) 
real  ch 
integer  k 

ch=h/4 . 

do  221  k= 1 , neqn 
f 5 ( k ) ~y ( k ) +ch*yp( k ) 
call  f(t+ch,f5,fl) 

ch=3. *h/32. 
do  222  k= 1 , neqn 
f5(k)=y(k)+ch*(yp(k)+3.*fl(k) ) 
call  f ( t+3 . *h/8 . , f 5 , f 2 ) 

ch=h/219?  . 
do  223  k= 1 , neqn 

f5(k)=y(k)+ch*( 1932. *yp(k)+( 7296. *f2(k) -7200. *fl(k) ) ) 
call  f(t+12.*h/13. , f 5  ,  f 3  ) 

ch  =  h/4  104 . 
do  224  k= 1 , neqn 

f 5 ( k ) =y ( k ) +ch* ((8341. *yp(k) -845. *f3(k))+( 29440. *f2(k) 
-32832. *f 1  (k) ) ) 
call  f ( t  +  h , f 5 , f 4 ) 

ch=h/20520 . 
do  225  k=l,neqn 

fl(k)=y(k)+ch*( ( -6080. *yp(k)  +  ( 9295. *  f 3 ( k ) 

-5643.*f4(k)))+(41040.*fl(k)-28352.*f2(k))) 
call  f ( t  +  h/2 .  , f  1  , f 5 ) 

ch=h/76 18050 . 
do  230  k= 1 , neqn 

3 (k ) =y <  k ) +ch* (  (902880. *yp(k)  +  (3855735.*f3(k) 

-1371249. *f4(k)))+(395 366 4. *f2(k)+277020.*f5(k)M 

return 

end 
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This  function  performs  a  table  look-up  to  determine 
the  rate  coefficient  at  a  given  energy,  en,  and  for  a 
given  process,  p. 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C************************************************************ 


& 


real  function  rate(p,en) 
real  en , on ( 27 ) , th ( 27 > , rates ( 27 , 
length , tai 1 ( 8 ) 
integer  p,n 
character*?  file(27) 
common/oncf  f /on , file 
common / f  req/th .rates, tail 
common/recom/y , temp, length 


561,2) , y ( 6 ) , temp , 


if ( ( on( p) .  eq  , 
rate=0 . 
return 
endi  f 


0. ) .or. (en.le.O. ) )then 


if( ( p. eq. 2 1 ) .or. (p.eq.23) 
rate: rates ( p , 1 , 2 ) 
return 


or. (p.eq.24! ) then 
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18 
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20 
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21 
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J 

23 

hdet 

$ 

24 

adet 
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e- 

25 
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26 

pdiss 

27 

V-T  relax 

'M 

JJ3 
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endi  f 


if{ (p.eq. 15) .or. (p.eq.16) )then 
i f ( en . gt . 1 . e7 ) then 
rate=0. 
else 

rate=8.75e-27/en**(4.5)*y( 1) 
endif 
return 
endif 


i f ( p . eq . 17 ) then 

rate=1.56e6*sqrt( temp/ 2 . ) /length 
return 
endif 


i f ( p . eq . 2  7 ) then 

rate=40./8.*1.5e-10*exp(-. 1603/ temp ) 
return 
endif 


i f ( en.lt. .01) then 
rate=0 . 
return 
endif 


i f ( en . le . 10 . ) then 

n= int ( en/ . 05+ . 5 ) 
goto  1 
endif 


if(en. le. 100. )then 

n=int(en/. 5+181 .5) 
goto  1 
end  i  f 


if(en.le.l000. )then 

n=int(en/5+361.5) 
i f ( n . eq . 56 1 )  n  =  560 
goto  1 
end  i  f 


rate=rates(p, 561,2) 
return 


rate=(en-rates(p,n, 1 ) )/(rates(p,n+l , l)-rates(p,n, 1 ) )  * 
(rates(p,n+l ,2)-rates(p,n,2) )+rates(p,n,2) 


re  turn 
end 


This  subroutine  calculates  the  effect  of  wall 
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c  collisions  on  the  average  energy  of  the  source  by 

c  modifying  the  Thomson  distribution  by  a  survival 

c  probability.  It  also  calculates  the  source  strength, 

c 

c ************************************************************ 

subroutine  thomson ( ens , vol t , so ) 

real  ens, volt, tol,integral(3) , tl ( 200 ) ,a3(200) ,h(200) , 
&  fa(200) ,fc( 200) ,fb(200) ,s(200) ,1(200) ,fd,fe, si, 

&  s2 , v ( 8 ) , so , eps , epspl 


eps  is  machine  epsilon,  or  the  limit  of  machine  accuracy. 
It  is  calculated  to  be  used  as  a  minimum  allowed  error 
tolerance . 

tol =  1  .  e-2 
eps  =  1 . 

>  eps=.5*eps 

epspl =eps  + 1 . 
i f ( epspl . gt . 1 .) goto  10 
eps=20 . *eps 

do  100  j  =  1 , 3 

integral ( j ) =0 . 

do  1  i  =  1  , 200 
tl ( i ) =0. 
a3 ( i ) =0 . 
h ( i ) =0 . 
f  a ( i ) =0 . 
fc ( i ) =0 . 
f b ( i ) =0 . 
s ( i ) =0 . 

1( i )=0. 

con  t i nue 

fd  =  0  . 
fe  =  0  . 
sl=0. 
s  2  =  0  . 

do  2  i  =  1  ,  8 
v ( i ) =0 . 

i  =  1 

1 1  (  i ) -  1 0 . *  to  1 
a3 ( i ) =0 . 
h(i)= volt/2, 
f  a ( i ) =  f  t ( 0 .  ,j) 
fc ( i ) =f t ( h ( i ) , j ) 
f b ( i ) =f t ( volt , j ) 

s( i)=h(i)*(fa(i)+4.*fc(i)+fb(i) )/3. 

1  (  i  >  =  1  . 


i 


3  if ( i . le . 0 ) goto  100 

if ( tl ( i ) .lt.eps)tl(i)=eps 

fd=ft(a3(i)+h(i)/2. ,j) 
fe=ft(a3(i)+3.*h(i)/2. , j) 
al=h(i)*(fa(i)+4.*fd+fc(i) )/6. 
s2=h(i)*(fc(i)+4.*fe+fb(i> )/6. 
v ( 1 ) =a3 ( i ) 
v( 2  )  =fa( i ) 
v ( 3 ) =fc ( i ) 
v( 4 ) =fb( i ) 
v( 5 )=h( i  ) 
v ( 6 ) =tl ( i ) 
v( 7 ) =s{ i  ) 
v( 8 )  =  1 ( i  ) 


if ( abs ( s  l+s2-v ( 7 ) ) . It . v ( 6 ) ) then 

integral ( j ) = integral ( j ) +sl+s2 
else 

if ( v( 8 ) . ge . 200 ) then 

write ( 6 ,*)’ Integration  failed  in  Thomson’ 
stop 
else 
i  =  i  +  l 

a3(i)=v( 1 ) +v( 5} 
f a ( i ) =v ( 3 ) 
fc(i)=fe 
fb( i ) =v( 4 ) 
h(i)=v(5)/2. 
tl( i )=v(6)/2. 
s { i ) =s2 
1 ( i ) =v ( 8 ) + 1 . 
i  =  i  +  1 

a3 ( i ) = v ( 1  ) 
fa( i ) =  v( 2 ) 
f c  (  i  )  =  fd 
fb( i ) =v( 3 ) 
h(i)=h(i— 1) 

1 1 ( i ) =  tl ( i-1 ) 
s ( i ) =s 1 
1  < i  )  =1  < i-  1 ) 
endi  f 
endi  f 

goto  3 

100  continue 

ens  =  integral (  1  )/ integral ( 2 ) 
so= integral ( 2 )/ integral ( 3) 
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return 

end 


c  This  function  is  the  integrand  used  in  the  routine  thomson 
c 

real  function  ft(e,j) 
real  e,p,x 

i f ( j . eq . 3 ) then 

ft=l./(e+15.427)**2. 
goto  2 
endi  f 

x=fwall ( e ) 
y=fcoll ( e ) 
i f ( x . eq . 0 . ) then 
p=l  . 
else 

p=l  .  -x/ ( x  +  y ) 
endi  f 

ft=p/(e+15.427)**2. 
if ( j .eq. 1 ) ft=ft*e 

2  return 

end 


This  subroutine  outputs  various  pieces  of  information 
at  the  end  of  the  run  in  order  to  better  understand 
what  is  going  on  in  the  discharge. 


subroutine  endout(y) 

real  y(6) , telec , r ( 27 ) ,cp,cn, pin, pout, th(27) ,ens,so,no, 
&  nh,ml,m2,rates(27,561,2),telec,relectron,tail(8), 

&  distrib, area, vol,vp, volt, curr,ftail(0: 100) 

common  so,no,nh 
common /ener/ens , m 1 , m2 
common/ f  req/th, rates, tail 
common/ geom/ area , vol , vp 
common /power/ vo 1 t , cur r 
common /d is t/f tai 1 


telec =2 . *y ( 5 ) /3 . 

write ( 20, 1 ) y  < 1 ) ,y( 2) ,y( 3) ,y( 4 ) ,y(6) , telec 
format < //20x , ’ SUMMARY ’//' Electron  density 
*  ’H-  density 

&  ’H2+  density 


’ , IpelO. 3/ 
’ , IpelO. 3/ 
’ , IpelO. 3/ 
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'H  +  density  =  ’ 

*H2(v)  density  =  ’ 

'Electron  Temperature  = 

write ( 20 . 13 ) y ( 2 ) /y ( 1 ) 
format (/’ N-/Ne  ratio  =  IpelO. 3) 
wri te ( 20 , 12) no, nh ,so,ens 
format(/’H2  density  =  ’,lpel0.3/ 

'H  density  =  ’,lpel0.3/ 

’Source  strength  =  IpelO. 3/ 

'Average  energy  of  source  =  ’ ,0pf6.3) 


IpelO . 3/ 
IpelO . 3/ 
, Opf 6 . 3 ) 


Charged  particle  balance 


cp=y ( 3 ) +y84r 
cn=y ( 1 ) +y( 2 ) 
write(20,2)cp,cn 

format (//’CHARGED  PARTICLE  BALANCE'/ 
'Total  POSITIVE  charge  density 
'Total  NEGATIVE  charge  density 


’ , IpelO. 3/ 

’ , IpelO. 3/) 


Rate  coefficients 


’ ,r< 15) 
’ ,r(  16) 


’ ,r( 18) 
' ,r(20) 


write (20,4) 

format (// 'RATE  COEFFICIENTS  IN  STEADY  STATE’) 
do  3  i= 1 , 27 
r ( i ) =rate( i ,y( 5 ) ) 

wr ite ( 20 , 9 )’ Dissociative  Recombination  =  ’,r(12) 
write ( 20 , 9 )’ Vibration  =  ’,r(14) 

write ( 20 , 9 )’ Recombination  -  H2  =  ’,r(15) 

write ( 20 , 9 )’ Recombination  -  H  =  ’,r(16) 

wri te ( 20 , 9 )’ Dissociative  Attachment  =  ’ , r (  1 8 ) 

wri te ( 20 , 9 )' Momentum  Transfer  =  ’,r(20) 

write( 20 , 9 ) 'Mutual  Neutralization  =  ’,r(21) 

wri te ( 20 , 9 )’ Detachment  by  e  =  ’ , r ( 2 2 ) 

wr i te ( 20 , 9 )’ Detachment  by  H  =  ’,r(23) 

wri te ( 20 , 9 )’ Associative  Detachment  =  ’,r(24) 

wr i te ( 20 , 9 ) ’ Polar i zed  Dissociation  =  ’,r(26) 

wri te ( 20 , 9 ) ’ V-T  relaxation  =  ’,r(27) 

rel ec tron= ( tai 1 ( 1 ) + tai 1 ( 2 ) +tai 1 ( 3 ) +tai 1 ( 4 ) 

+  tail(5)+tail(6)+tail(7)+tail(8)  )  / y (  1) 
wr  i  te ( 20 , 9 )  ’ Total  electronic  -  tail  =  ’.relectron 

relectron=r( 1 )+r(2)+r(3)+r(4)+r(5)+r(6)+r(7)+r(8) 

+  r(  10  )  +r (  1 1 ) 

wr i te ( 20 , 9 ) ’ Total  electronic  -  MB  =  ' , relectron 

formatfa, IpelO. 3) 


’  ,  r  (  2  1  ) 
’  , r ( 22  ) 


’  ,  r  (  2  3  ) 
’  ,  r  (  2  4  ) 


, r ( 26  ) 
,  r  (  2  7  ) 


Power  balance 


write(20,5) 

format (//’ POWER  PARTITION’) 
pout=0 . 

p i n  =  ens  *  so/y (  1  ) 
x  =  0  . 

do  6  i = 1 , 8 
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c 

c 

c 

c 


x=x+r( i ) *th( i  ) 

x=(x+r( 10)*th( 10)+r(ll)*th( 11) )*no 
write ( 20 , 10 )’ e-beam  power  input 
format ( a , lpe 10 . 3 ) 

write ( 20 , 1 1 )’ power  into  electronic 

pout=pout+x 

x=no*th< 14 ) *r ( 14) 

write ( 20 , 1 1 )’ power  into  vibrational 

pout=pout+x 

x  =  y( 3 ) *y( 5 ) *r( 12 ) 

write ( 20 , 1 1 )’ power  into  diss  rec 

pout=pout+x 

x=no*th( 9 )  *r ( 9 ) 

write ( 20 , 1 1 )’ power  into  H2  ionization 

pout=pout+x 

x=nh*th( 13)*r( 13) 

wr ite ( 20 , 1 1 )’ power  into  H  ionization 

pout=pout+x 

x=y<3)*y(5)*r<15) 

write ( 20 , 1 1 )’ power  into  recombination 

pout=pout+x 

x=y( 4 ) *y( 5 ) *r ( 16 ) 

write ( 20 , 1 1 ) 1  power  into  recomb  of  H 

pout=pout+x 

x=y(6)*y(5)*r( 18) 

write ( 20 , 1 1 )' power  into  diss  att 

pout=pout+x 

x=2.*ml/m2*y{5) *no*r ( 20 ) 

write ( 20 , 1 1 )’ power  into  momentum  trans 

pout=pout+x 

x=y(2)*y(5)*r(22) 

write( 20 , 1 1 )’ power  into  detach  by  e 
pout=pout+x 

x=y(2)*nh/y( 1 ) *y ( 5 ) * r ( 23 ) 

write ( 20 , 1 1 )’ power  into  detach  by  H 

pout=pout+x 

x=y(2)*nh/y( 1 ) *y ( 5 ) * r ( 24 ) 

write ( 20 , 1 1 )’ power  into  assoc  det 

pout  =  pout.  +  x 

format ( a , lpe 10 . 3 , 2x , Opf 7 . 2 , ’%’ ) 


write(20,8)pin, pout 

f ormat (/’ Total  power  input  =’,lpel0.3/ 
’Total  power  output  =’,lpel0.3) 


( pm 


this  section  of  code  outputs  the  EEDF  into  the  file 
dist.out  in  the  correct  format  to  be  used  by  'graph’ 


*0 


open (unit=15,file=’dist.out’ ,  status =' new ’ ) 
write( 15,*) ’Electron  distribution ’ 
write(15,*)’ Energy ’ 
wr i te ( 1 5  ,  *  )  ’  ( eV ) ’ 
wri te ( 15 , * ) ’ f tot ( e ) ’ 


Page  102 


8 


’ , x , x/pin* 100 

’ , x , x/pin* 100 

& 

’ , x , x/pin* 100 

i 

v: 

v; 

’ , x , x/pin* 100 

■>, 

’ , x , x/pin* 100 

f.i 

ft 

’ , x , x/pin* 100 

& 

& 

’ , x , x/pin* 100 

i 

WO 

’ , x , x/pin* 100 

■ 

’ , x , x/pin* 100 

3 

£< 

’ , x , x/pin* 100 

% 

W 

4 

'  , x , x/pin* 1 00 

•.y 

8 

’  , x , x/pin* 100 

8 

1 ) 

qrt ( 4/3 . 14159) *y ( 1 ) /telec** ( 3  . /2  .  ) 
exp( -en/telec ) 
istrib+f tail ( i ) 

* ) en ,dist 
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